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Development and correctness analysis of the mathematical model of transport and suspension 
sedimentation depending on bottom relief variation ° 


A. I. Sukhinov’, V. V. Sidoryakina’ 


' Don State Technical University, Rostov-on-Don, Russian Federation 
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Ilocrpoenne u ucc1eqOBaHhe KOPpeKTHOCTH MaTeMATHUECKON MOeJIM TpancnopTa HW OcaxKeHHnA B3Becei 
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* Taranporcknii nuctutyt umeun A. II. Uexoza (cpunnam) PrOY (PHHX), r. Taranpor, Poccuiickas Megzepannsa 


Introduction. The paper is devoted to the study on the three- 
dimensional model of transport and suspension sedimentation 
in the coastal area due to changes in the bottom relief. The 
model considers the following processes: advective transfer 
caused by the aquatic medium motion, micro-turbulent diffu- 
sion, and gravity sedimentation of suspended particles, as well 
as the bottom geometry variation caused by the particle set- 
tling or bottom sediment rising. The work objective was to 
conduct an analytical study of the correctness of the initial- 
boundary value problem corresponding to the constructed 
model. 

Materials and Methods. The change in the bottom relief aids 
in solution to the initial-boundary value problem for a parabol- 
ic equation with the lowest derivatives in a domain whose 
geometry depends on the desired function of the solution, 
which in general leads to a nonlinear formulation of the prob- 
lem. The model is linearized on the time grid due to the “freez- 
ing” of the bottom relief within a single step in time and the 
subsequent recalculation of the bottom surface function on the 
basis of the changed function of the suspension concentration, 
as well as a possible change in the velocity vector of the aquat- 
ic medium. 

Research Results. For the linearized problem, a quadratic 
functional is constructed, and the uniqueness of the solution to 
the corresponding initial boundary value problem is proved 
within the limits of an unspecified time step. On the basis of 
the quadratic functional transformation, we obtain a prior es- 
timate of the solution norm in the functional space L2 as a 
function of the integral time estimates of the right side, and the 


initial condition. Thus, the stability of the solution to the initial 


Beedenue. Hactostlaa padota MocBalleHa MCCIeOBaHHIO 
IIPOCTpaHCTBCHHO-TpexMepHOH MOesIM TpaHcioptTa MU oca- 
%KCHHA B3BECH B MpHOpexHOM 30HeE C yYeTOM H3MeHeHHA 
pembeda gHa. Moyemb y4uTbIBaeT cIeytyroulve Mporeccsl: 
aJIBEKTUBHbIM TepeHoc, OOyCJIOBJICHHbIM JIBUXKeEHHeM BOJHOM 
Cpelbl, MUKpOTypOyIeHTHy!o Wuddy3ulo U rpaBuTallMOHHoe 
OcaxkK]IeHve 4aCTHI, B3BECH, a TakxKe M3MeHeHvWe TeOMeTpHU 
J{Ha, BbI3BAaHHOe OCaKeHHeM 4aCTHI B3BeECH WIM MObeEMOM 
YacTUI JJOHHbIX OTIIOXKeCHHH. 

Llenbio pabomol Aéi1AIOCb TpoBeeHve aHaMTU4eCKOTO UC- 
Cyle]OBaHHA KOPPeKTHOCTH HavasIbHO-KpaeBou 3aja4dH, COOT- 
BeETCTBY}FOINeM MOCTpOeHHOM MOEN. 

Mamepuaaoi u Memoooi. U3mMeHeHne pesbeda Ha IpMBOMT K 
HeOOXOXMMOCTH pelllaTbh HavasIbHO-KpaeByloO 3aya4dy lA 
ypaBHeHua Mapadommueckoro Tula C MJIayWIMMM Mpon3Bo]- 
HbIMH B OOIACTH, TEOMeETPHA KOTOPOM 3aBHCHT OT HCKOMOM 
(YHKUMU pellieHua, YTO MPUBOMT, B OOMIeEM Ciyyae, K HeJIN- 
HeMHOM MOCTaHOBKe 3aa4nH. BeinosHeHa JMHeapu3allua MO- 
jjeIM Ha BPCMeHHOM CeTKe 3a CYeT «3AMOPAXKUBAHHA» PelIbe- 
(ba WHa B ipeyjeswax OJHOrO Iara 10 BpeMeHH HU Wocuesyrolle- 
TO Mepecueta (PyHKIMH MOBepXHOCTH Ha Ha OCHOBe H3Me- 
HUBIIeHcaA (YHKUMH KOHI[CHTpallM B3BeINICHHOTO BeLIeCTBA, 
a TaKXKe BO3MO2KHOTO H3MeCHeEHHA BeKTOPa CKOPOCTH JIBYDKe- 
HHA BOAHOU CpeJBl. 

Pe3zyibmamoi ucciedoeanus. JIA IMHeapv30BaHHOu 3aqa4qn 
TOCTpoOeH KBal{paTHyHbI MyHKUMOHAaI WU 9HepreTHYeCKHM 
MeTOJIOM J{OKa3aHa eCJMHCTBCHHOCTb pellleHHA COOTBETCTBY- 
role HavaIbHO-KpaeBOU 3aa4u B IIpeyjeiax MpOU3BOJIBHOrO 
Iara 10 BpeMeHH. Ha ocHoBe TpeoOpa30BaHHA KBapaTH4HO- 
ro PyHKWMOHaIa MOWy4eHa allpuopHasd OLleHKa HOPMBI peliie- 
HUA B PYHKIMMOHAIbHOM WpoctpaHctTBe Ly B 3aBHCHMOCTH OT 
WHTerpasIbHbIX OIeCHOK MO BpeMeHH IpaBow YacTu, rpaHu4- 
HbIX YCJIOBHM WM HavasIbHOrO YyCJIOBHA, HU, TaKHM OOpa30M, 
JJOKa3aHa YCTOMUMBOCTh pellieHHaA HCXOAHOM 3aa4n pu u3- 
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problem from the change of the initial and boundary condi- 
tions, the right-hand side function, is established. 

Discussion and Conclusions. The model can be of value for 
predicting the spread of contaminants and changes in the bot- 
tom topography, both under an anthropogenic impact and due 
to the natural processes in the coastal area. 


Keywords: coastal systems, mathematical model, diffusion- 
convection problems of suspension sedimentation, bottom 
relief change, uniqueness of solution, and stability of initial- 
boundary value problem. 


For citation: A.J. Sukhinov, V.V. Sidoryakina.. Development 
and correctness analysis of the mathematical model of 
transport and suspension sedimentation depending on bottom 
relief variation. Vestnik of DSTU, 2018, vol. 18, no. 4, 
pp.350-361. https://doi.org/10.23947/1992-5980-2018-18-4- 
350-361 


MeHeHHH HavasIbHOTO WM rpaHW4HbIx ycOBHH, PyHKUMM Mpa- 
BOH 4acTH. 
O6cyacdenue u 3aknmouenua. Mojeltb MOxKeT TpeCTaBsATb 


I[CHHOCTh MPH MpOrHo3e paciipocTpaHeHuad 3arpa3HeHH U 
VW3MeHeHHA pebeda WHa, KaK IPH AHTPOMOreHHOM BO3JeH- 
CTBHU, TaK MW B CHJIy CCTECTBCHHO MpOTeKarollux MIpupOHbIx 


IIPOMeCccoB B IpHOpexHO 30He. 


KuroueBbie C10Ba: TpHOpexKHbIe CHCTeMbI, MAaTeMaTHYeCKaA 
MOJIeIb, 3a7a4u JHdPy3HH-KOHBeCKUMM OCaxKICHUA B3BeLIeH- 
HOTO BelllecTBa, H3MeHeHHe pelbeda Ha, CAMHCTBCHHOCTb 
pelleHua UW yCTOMYMBOCTh HavasIbHO-KpaeBoHu 3aa4H. 


Oopazeu ona UUumuposeanua: CyxuHos, A. VU. locrpoennue u 
ucciieqOBaHve KOPpeKTHOCTH MaTeMaTHYecKOH MOseuM 
TpaHcnopTa H OcaxKeHHA B3BeCeH C YUCTOM H3MeCHeHHA pelib- 
eda qua / A. UW. Cyxunos, B. B. Cugopsxknua // BectHuk Jon. 
roc. TexH. yH-Ta. — 2018. — T. 18, Ne 4. — C.350-361. 
https://doi.org/10.23947/1992-5980-2018-18-4-350-361 


Introduction. The aquatic habitat protection [1—2] is one of the most important factors that determine the inte- 
grated research development of the coastal areas. Damage control over the natural processes, such as pollution, sedi- 
mentation, and depletion of water areas, leads to necessity for studying all aspects that affect changes in coastal waters. 
Maintenance of water bodies in proper condition and timely intervention in its operation mode is directly related to the 
increase in port capacity and the efficient development of the coastal infrastructure (ensuring an accessway to the berths 
of ships with a low landing; desilting and aquatic vegetation clearing of the coastal strip; etc.) [3-5]. As a rule, research 
practice in this field requires the construction of mathematical models that are as close as possible to real processes 
[6-11]. 

A continuous mathematical model describing spatial-three-dimensional processes associated with transport and 
gravitational suspension sedimentation in the aquatic medium with varying bottom relief 1s considered. This model 
takes into account micro-turbulent diffusion and advective transfer of suspensions, the effect of gravity on suspension, 





the presence of the bottom and a free surface, and a bottom contour variation. 

The suspension transport model enables to study the hydrophysical processes of aquatic systems, to predict the 
dynamics of the bottom surface change based on the description of the lifting, transport, sedimentation, changes in the 
concentration of suspension [12—13]. The uniqueness of the solution to the corresponding initial-boundary value prob- 
lem is proved, and a prior estimate of the solution norm is obtained depending on the integral estimation of the right- 
hand side, boundary conditions, and the initial condition. 

Materials and Methods. Continuous 3D model of suspension diffusion-convection and the corresponding initial 
boundary value problem. Consider a continuous mathematical model of sediment spreading in the aqueous media 
considering diffusion and convection of suspension, gravity action on suspension, presence of the bottom and a free 
surface. We will use Oxyz Cartesian coordinate system where Ox axis passes along the nonperturbed water surface 


and is directed toward the sea, and Oz axis is directed vertically downwards. Assume that h=H-+n is the total water 


depth, m; H is depth with undisturbed water surface, m; n is elevation of the free surface relative to the geoid (sea lev- 
el), m (Fig. 1). 








~ Geoid — 4; Freee surface 
* 2 
if x 
, > 
: y . > 
jz 
Bottom __ _ 
——— ey ee ee 
z = H(x,y) 


Fig. 1. Introduction of Oxyz coordinate system Oxyz 
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Suppose that in G = {0 S15b,0sysl,05257 (x, y)} closure region, there are suspensions which have 
eC, ys z,t) concentration at (x, y,z) point and at ¢ time, mg/l; ¢ is temporary variable, sec. We will also use 
L= so DAX H (x,y) notation. 


The behavior of the suspended particles will be described by the following system of equations: 


de, Owe) , (ve) | O((w male) ES a ( +r, 


a a | by az a Oy | a & (1) 
OH € 

—=-—w,c, 

Ot po * 


where u, v, w are components of U fluid velocity, m/s; w, is hydraulic size or sedimentation rate, m/s; UW, 1, are coef- 
ficients of the horizontal and vertical turbulent diffusion of particles, respectively, m’/s; F is power of particle sources; € 
is porosity of bottom materials. 

Summands on the left side (except for the time derivative) of the first equation of the system (1) describe the 
advective particle transport due to the inertial motion of the aqueous media, as well as sedimentation under the action of 
gravity. The summands on the right side describe the suspension diffusion. The vertical diffusion coefficient 1s chosen 
different from the horizontal diffusion coefficient due to the fact that the effect of difference between these coefficients 
is often observed in various media and can be caused by various factors. 

As G region, we consider ABCDAOC,D, “parallelepiped” “skewed” to the shore, whose 4 OC,D, upper base 


lies on (z = 0) free surface, and (z = A(x, y)) part of the bottom surface is its lower base. Suppose S is G surface, /i 


is the outward normal to the surface of the “skewed parallelepiped”. We assume the given U’ as the fluid velocity on 
G side surfaces. Complete with the boundary conditions of first kind for the particle concentration function, this allows 
determining the suspension flow both towards the coast and along the coast (Fig. 2). 


Cc; Coast 














Free surface 






Particle flow 
precipitates 









Particle flow 


precipitates 
D 





Bottom 


Fig. 2. Solution area for suspension transport 


Add the initial and boundary conditions (assuming that the sedimentation is irreversible) to the system (1) 
As the initial conditions at ¢ = 0 time, we accept 


c(x, y, z,0) =; (x, y,z) ; (2) 
H(x,y,0) =H, (x,y). (3) 
We set boundary conditions on ABCDAOC,D, faces (we set suspended flows both towards the coast and 
along the coast): 
- on the faces S,=AAOB (x=0, OS y<L, 0<z<L,) , S,=AADD (y=L,, 0Sx<L, 0<z<L,) and 
S,=BOCC (y=0, 0<x<L,, 0<z<L) 
c=c ,Tme c =Cc Ket) t<[0,7]; (4) 
.) and S,=AOCD, (z=0, 0<x<L, 0<y<L,) 
C=O: (5) 
0<y<L,) 


- on the faces S,=DD,C,C (x=L,, 0S y<L,, O<z<L 


- on the surface S, = ABCD (z =H(x,y,t), OS x<L 


x? 
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w w 
te, WIM nl ae (6) 
On LL, Oz U, 


The boundary condition (5) occurs with a relatively small slope of the bottom: 


| oH Y 
max ,/|/ ——|] +| —| <«l. 
So Ox Oy 


The following condition of the solution domain nondegeneracy is set up for all (x, y,t) at which the initial 


boundary value problem is formulated: 

H (x, y,t) 2h, = const > 0, O<t<T. (7) 
When studying combined models of sediment and suspension transport, it is possible to increase the concentra- 
tion of suspended particles in the bottom layer due to the rising bottom sediment particles if the shear stress exceeds of a 
certain critical value is exceeded[13—16]. Then, instead of the boundary condition (6), we will consider the boundary 

condition of the following form 
Oc _ 
o 
Linearization of the initial-boundary value problem of transport and suspension sedimentation. To create a line- 


ac, a&=const > 0. (8) 


arized model on 0<t<T time interval, we construct a uniform grid wt with a step t, that is, a set of points 
oO. =40 =nt, N= Oli JV, Near. 





Fig. 3.Construction of time grid 


oc”) (x, y,z,¢,,)and H (") (x, y,t,, ) functions are determined at each step of w, time grid. If n=1, then func- 
tions of the initial condition will suffice to cl” (%95251,), df ) (x,y), viz cl (x, y,z,0)=c, (x,y,z) , 
H" (x, y,t,) =H, (x,y) respectively. But ifn = 2,...,N , then oc) (x, y,z,t,,) = cv) (x, y,z,t,,) functions are assumed 
to be known, since the problem (1)—(6) for the previous t_, <¢<t,, time interval is supposed to be solved. 


We write the system (1) on t., <f¢<t, interval in the form: 


ate) so) ltrna ie) fh oe) af, a, 














at a. o az a? ay | al” (0) 
eH”) & (n) 
=-—Ww.c 
Ot po ° 
and complete it with the initial conditions: 
cl) (x, V,z, t, ) =, (x, y,z),c” (x, V2. is) =) (x, V2 ), = 2 ccc: (10) 
H") (57, = H, (x,y), H” (x, y.t,-1 ) = Fer) (x, y.t,-1 ), n=2.,..., N. (11) 


The boundary conditions (4)—(6) are assumed to be fulfilled for all ¢,, <¢<¢, time intervals. 
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By defining c")(x,y,z,t,,)=cl""(x,y,z,t,,) function on t,,<¢<t, time interval, we can find 
H” (x, Vis) function. For this end, we integrate both members of the second equation of the system (9) over 


t_, <t<t, variable. We will get 





t t 
n yy”) n 
| : dt = ~~ wy, | ol" dt. (12) 
tn-1 ot p ty 
From the equality (12), it is not difficult to get 
N_ i'n 
H” =H’) ices | ol" dt. (13) 
n=l ¢ ' 


We introduce G,,= {0 <x<L,0<y<L,0<z<H” (x, y,t,, )} domain at each ¢,, <t <t, time step. 


We have a chain of linear initial-boundary value problems for each time layer, where the system of the type 

















acl”) a(uc"” a(vel”) a((w+ w, Je”) ac atl) a acl”) (14) 
+————- + + ——_—_————=, | +F, 
Ot Ox Oy Oz Ox Oy Oz Oz 
(xy,z)eG_,, G, — O<y<L,, 0<z<H""(x,y,t,.)}; 
(n) _ HY) _& une 
H med | ea | Po eee (15) 
is considered for ¢,, <¢<¢, interval with the initial conditions: 
o”) (x,y, z,t,,)=c"" eae, Z,t,1)s (16) 
HH") (x, Vet 4 ) =e (x, Vit 4 ) (17) 


Note that at each time step, the boundary surfaces will change (except S, face). Considering ¢t,_, <t<¢, time 


interval, we set the boundary conditions on the edges of G_, domain: 


i(¥=0, O<y<Z,, O0<z<H"(0,y,t,)) , S,,.(y=L,, OSxSL,, OS z<H(x,L,,t,,)) and 


- on 8, er 


S,.i(y=0, OSx<L,, 0Sz<H""(x,0,t,,)) faces 
co”) =e" ,rne c’ =e" (x,y,z,t), te[t,,.t,]; (18) 
-on S,,,(x=L,, OS ySL,, OSz<H""(L,,y,t,,)) u S,,,(2=0, OSxSL,, OSySL,)=AOCD, faces 
ce” =0; (19) 
-on S,,,(z=H""(x,y,t,,), OS XSL, 0S y<L,) surface 


Aee) w, (n) ac”) Ww, (x) 











=—-—c'’ wm =—-—C (20) 
on LL, Oz LL, 
The boundary condition (8) will be replaced by the following 
(n) 
acl”), a=const > 0. (21) 
OZ 


Thus, it is supposed that the bottom relief within this time step, when calculating the distribution of suspension 


concentrations, does not change and is taken from the previous time layer. First of all, at this ¢, <¢<t, time step, the 


initial-boundary value problem for the convection-diffusion equation (14) with H'” fixed bottom relief function is 


solved, and only then the update (recomputation) of H”) relief function is performed in accordance with the equality 
(15). 

The determination of the conditions of existence, uniqueness and continuous dependence of the solution on the 
input problem data is carried out on a fixed time layer under these assumptions and subject to the condition (7). 

The authors do not plan to study the existence of solutions to the initial-boundary value problems (14) - (20) 
and (14) - (19), (21) in this paper. Questions of the existence of solutions to the initial-boundary value problems for 
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parabolic equations with lower derivatives (diffusion-convection equations) are considered, for example, in the mono- 
graphs [17-18]. 

Research Results. Investigating uniqueness of the solution to the initial-boundary problem of suspension 
transport. 
Consider the initial boundary value problem (14) - (20) formulated for the arbitrary t,,<¢<t, time layer. 


Multiply the left and right member of equation (14) by c’” function and get: 


oo i A(uc'” | re a(ve'”') 4 a((w+w, Je”) _ uc”) (a in Oc”) |. ol”) O Ant |. ol F (22) 











Ot Ox Oy Oz Ox? éy OZ 


The left member of the equality (22) can be transformed as follows: 


122 uf ee") ae”) alwende)) 12(eF sag (uz. 














0 2 ot 
Ot Ox Oy . Zz (23) 
a(cl”’ b 
rele") +5aiv((c) C 
2 Ot 2 
where U =|lu, v, w+ w, || - 
With regard to (23), the equation (22) will be written as 
2 

a(c” _ a(n) a2 A() (n) 
. | | +5aiv((c) U)=n,c" aie — a a +00) = a 1, —— a +0 F. (24) 
2 ot 2 Ox Oy Oz Oz 


Then we integrate both members of the equation (24) over ¢_, <¢<t, interval, and, after that, over the spatial 


variables in G_, domain. In the first term, the order of integration is changed due to the Fubini theorem [19]. We obtain 
2 
tr O ao ‘ 1 os 
eel Gal J S| [Nae((e?) Gag, far 
G,-1 2 a Ot 4 2 G,_1 


lw ol | fee] (fe ate “Jac, i (25) 
j[ferra. 


The first term on the left side of the equation (25) is obviously equal to 


» aly 
III; ple, a6, = [JJ 5((e) (s9.204,)-(e) (7.28.0) dG, (26) 





Next, we turn to the transformation of the second term of the left-hand side of the equality (25). Considering 
the Gauss-Ostrogradsky formula and the boundary conditions (18) - (20), it can be written as [20]: 


ic i aiv((c" oa, je=F 4 (9 (6) (Flare or Ff if Cw dy a 


-1 


*| ff (ey (Tadeo Ff I ( yo Alaa} F[ 2) ud (27) 


vee Ae 


-Fi( se y nde Jo 5 fl (C) vd +5 ff (e”Y wavy a 


ta \S3n4 


where U" is the known velocity of the aquatic medium on the faces where the boundary conditions of the first kind 


and S,_, top cover on which the suspension concentra- 


-l 


are specified; in fact, these are all side faces, except for S, , > 


tion is zero, and therefore the flows through them are zero. 
Let us turn to the transformation of the right side of the equation (25). The following equality occurs 
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Ae) Ao” a Ae”? 
iM" TH, _— HM, - dG, = 

ax. Oy\ oy OZ OZ 

ol”) () () 

= al ree +L, a c") is mc Ou, ue dG..- (28) 

Ox Ox oy Oy OZ OZ 
acl”) ac!” ac”) 

— ——| + + dG... 

‘ Mtl ay J ay JO a a 


Ac”) os ac” fn dc” 


Suppose O=!0,0, =suc) —, wc" Be 
ppose O={0.,0,,0.} " ape a 
































} Then, in virtue of the Gauss-Ostrogradsky 


theorem, we have: 











a) (n) 
I) LL, 2am e Jon 2{e" = }oZ fer, - “es. ={I] diva - 


= J O, dxdz + J O dydz+ I O, dxdz + I O dydz + I O. dxdy + ‘Tf qa (29) 
. - iI O, dxdz + ai O, dxdz + ii O.dydz+ I cc ae 


Transforming each term fom the aphenned side of 09) subject to tic conditions on the boundary (18) - (20), 


we obtain 
ol”) (”) (n) 
fy) 2{er% Jou 2 fem : }S{ems : Jeo. 
x x Oy y Z z (30) 


=I! C 1, Lede {| C 1, Sede J C 4, Save ie Le OV dxdy. 


S31 —11 











In virtue of 06), (28), (29) and (30), the equality (25) aes on form 


Le) (2 ui Hewes, Sat 
i(abe V v+e'y, Jane beo fl (He) | pe oe 
i otto OE) (5) FE] ey 

Efile) orate. jf ffferrac. fs 


The identity (31) will be fundamental under studying the uniqueness and obtaining a prior estimate of the solu- 
tion norm of the initial boundary value problem (14) - (20). In case of replacing the boundary condition (20) with the 
boundary condition (21), the quadratic functional (31) changes as follows: 


zilery vom lg Hey ween, Soa fa 
ater pi, ey lyGe “ys yen, te 
a) 1 v1 } |. dt = 


fle NV (ys20ty- “ael(G FAG, 


tnt Gh 


(31) 











(32) 











Suppose that the equation (14) with the same conditions (16) - (20) satisfy two different solutions to 
C=C (x, ye. z,t), C=C (x, ys z,t) problem. For their C= c, —c, difference, the following initial-boundary problem is 


valid: 
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ae, (ue) ave) | Al(we mje) _ (Se ( ) (33) 


SS + — ——' |, 
a. oe az "la ay | al?’ & 
c(x,,z,0) = 0), (x,y,z) eGs : (34) 
- on Si 1? ae > ee > Dig > Da faces 
€=c -c =0; (35) 
-on S,,, surface 
Ww Ww 
OH, H, 


For c function, the equality (33) will take the form considering the equalities (34)—(36) 


fle x YZ, )dG,. +> S| ff mea |e 


-l1 S6 1 


[| ‘). n(Z) nS) Js |a-0 


Since w, >0 and other known values under the sign of integrals are positive uh> 0, uv> 0, then the equality 


(37) 


(36) 1s satisfied only under the condition 
E(x, y, z,t) = 0), (x, y,Z) eG.4, £45051, (38) 
which completes the proof of the uniqueness of the initial-boundary value problem (14) - (20) solution. 


In case of replacing the boundary condition (20) by the relation (21), instead of the expression (37), we obtain 
the following equality 


Seat 


2 [flO (n2t,)d6 aif (+ sa, Jd | 


-] (39) 
Oc Oc Oc 
+f] | (Zo (Z) on (SJ fe, aro 
We require the fulfillment of the inequality 

l 

ae —ap, 20, (x, y,Z) SS45 Cay st sh 
or 

Ww 
a<—, (x, yz) Site Gaus (40) 





2H, 
then all the terms in the equation (39) are nonnegative, and zero equality is possible if and only if 
C (x, y, z,t) = 0), (x yz z) EeG,, t_,<t<t,, that means the solution uniqueness and in this case. 

Reasoning is similarly repeated for all layers of w, time grid. The modification of the boundary conditions as- 
sociated with the continuous change in the bottom relief depending on the time variable requires additional study and is 
going beyond the scope of this article. 

Theorem. Suppose we are given a system of equations 








At ax Oz "| a? a’ | Oz 
(2 ” , G. g ateeech, 0<y<L,,0<z<H""(x,y,t,,)t, 
H” =H << Shera n =1,2,...,N 

p n=l 4 


inQ,,=G,,x(t,,<t<t,), G, = (0 <x<L,,0<y<L,,0<z< no (x, y,t,, )). simply connected domain with a 


sufficiently smooth boundary defined by the smoothness of z = H"”” (x, y) , O<x<L,, O<y<L, function with the 
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initial and boundary conditions (16) - (20). Let the functions of ee, y,z,t,,) solution, the velocity vector of 


r . . n—- .  e,e Te . * 
u,V,wt w, | aquatic medium, cv) (x, V2. t,,) initial condition, F (x, y, z,t) , right member of c (x, yz aa) bounda- 








ry condition, HW, = LL, (z), (x, y,Z) EG_, coefficient of u, =p, (z), (x, y,z) €G_, vertical turbulent exchange satisfy the 


following smoothness conditions: 











ol” (x, y,z,t,,)€C (Q,,)A Co. }, grad cc” € C(Q,,, ), U,V,w+w, | eC'(Q,,)0 Cio. ). 
ch (x, y,2t,,)€C(G,, ] 9 F(x, y,z,t)€C(Q,,) ? LU, (x, y,z)eC (G_,)C(G,,) e 
c’ (x, y,z,t)<¢C(S,,)x[z,, <¢s¢,], SS =G_\G, 
<ec((0<x<1, 0a 9 = 7. z=H"" (x,y))x[t,, <tSt, , as well as c'(x,y,z,0)=c,(x,y,z), 

n 
(x,y,z) ES,,\(0<x<L,0<y<L,, z= (ay); St CG, 


(0 <x<Lb,0<y<L, z=H Co (x, y)), conditions of consistency of the boundary and initial conditions, then the 


solution to this problem exists and is unique. 
Comment. In case of replacing the boundary condition (20) with the boundary condition (21), the inequality 
(40) should be added as a sufficient condition for the previous theorem. 

Studying the continuous dependence of the solutions to the initial-boundary value problem of suspension 
transport on the initial, boundary conditions and the right-hand side function. The next stage is connected with 
the study of the continuous solution dependence on the functions of the right-hand side, boundary and initial conditions 
for the system (14)-(15). 

Suppose that 
c >c, =const > 0, (41) 
02247, 02 ysl, 0S7<H” (a9t 4), by SPSe. 
For convenience, we introduce the notations: union of all parts of the lateral cylindrical surface (boundaries of 


G,,.; region) is denoted as S_,, and the lower base of region—as G_,— S,__, . In virtue of the smoothness conditions 


wn—l ? 


listed under the above theorem, extrema of functions on the bounded closed sets are reached: 


ac”) Ae”) 
) (42) 











) 
Sn 














Mg max {|e \, MA a vax 


Ox Oy 
v| e M,,-1 = min {Hy sb, \ 


9 








M = max { M = max {|u 
3,n-l Sent u,}. 4,n-1 Son1[ty1 StSt, | 


We will focus on the equation (31) if the boundary condition (20) is used, and on the equality (32) in case of 
the boundary condition (21). Evoking Friedrichs inequality, we have a chain of inequalities: 


acl") aol") \ acl") 
Cc Cc Cc 
—|+ + dG, ,2 
i f Ox sf ay [ Or n-l 
2 2 2 
acl”) ac” ac”) 
= minut} If Ay | ( ta | \eCet= (43) 


We turn to the equation (26) from which, in virtue of (42) and (43), we obtain the inequality: 
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feral el Zetz3* Gey fle 


2 
bn-l Git 


3 | {J Ww, iC ) iy dt < {fl G, dG. +My | {| (c" ) dS\_, h + (44) 


ty 
Chl S; nl bl 


t, " t, 
2M My, | iI elas, Jrvaa. | iielac,.. 2 
tha \ Sent tha \ Gai 


From inequality (44), there are two inequalities 


ite” dG, < [[fectG +My. | [fey as, jus 


th Den 




















(45) 
t, Ln 
OM, (Ma. I “|S. Joa | i[lelac,.. 2 
bet \ Seat bya \ Ga 
and 
f(e") dG, <M, ( [[Je.2dG,, + 
” 7 (46) 
+M,,, . {| (c’) 5, }-20,M,. . [fe is, }2a,. | ci jac} 
tra \ Sent tay \ Sea tt \ G4 
-1 
where M,,, = a 


= TU res 5 + ear 
2M ae L. L, ( HH" 


The inequalities obtained imply the continuous dependence (stability) of the solution to the problem (14) - (20) 
on the functions of the initial condition, the boundary conditions and the right-hand side, in Z, norm for any instant of 


0<7T <-+0 time, and also in ZL, time-integral norm. 


Obviously, if the inequality (45) and the theorem condition are satisfied, the initial-boundary problem (14) - 
(19), (20) will also have a solution that depends continuously on the functions of the initial condition, the boundary 
conditions and the right-hand side in the corresponding norms. 

Discussion and Conclusions. Novelty of the proposed non-stationary spatial-three-dimensional mathematical 
model of suspension transport lies in the fact that, alongside with considering the processes of advective transfer, micro- 
turbulent diffusion and gravity sedimentation of suspended particles, the model describes the change in bottom geome- 
try caused by the particle settling or bottom sediment rising. 

The linearization of the corresponding initial-boundary problem on the time grid is carried out, and the conditions 





for the uniqueness of the solution to the initial-boundary problem and continuous dependence on the input data — on the 
functions of the initial condition, boundary conditions, and the right-hand side in LZ, Hilbert space norm in ZL, time 


integral norm for two variants of boundary conditions are obtained for the arbitrary ¢,, <t<¢, time step. 
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Uncaennoe MOseMpoOBanHe Nomepe4Horo OOTeKaHHA MpOJIeTHBIX CTpOeHHii 6aIOUHBIX MOCTOB 


FO. A. Focrees’, A. JI. OGyxoscxnii’, C. J. Canenko>™” 


'? HopocHOupcKHi TOcyapCTBeHHBIM TeEXHH4eCKHi yHuBepcurer, r. HopocubupcK, Poccniicxad Dezepanna 


Introduction. The technique of numerical modeling of the 
transverse flow over span structures of bridges on the basis of 
the two-dimensional URANS (Unsteady Reynolds-averaged 
Navier-Stokes) approach used in the modern methods and 
software packages for computational fluid dynamics 1s verified. 
The work objective was debugging and experimental substanti- 
ation of this technique with the use of the database on the aero- 
dynamic characteristics of the cross-sections of span structures 
of girder bridges of standard shapes pre-developed by the au- 
thors. 

Materials and Methods. A numerical simulation of the trans- 
verse flow of low-turbulent (smooth) and turbulent air flows 
around the bridge structures in a range of practically interesting 
attack angles is carried out. SST k — w. turbulence model was 
used as the closing one. The technique was preliminarily tested 
on the check problem for the flow of the rectangular cross- 
section beams. Calculations were carried out using the licensed 
ANSYS software. 

Research Results. The calculated dependences on the attack 
angle of the aerodynamic coefficients of forces (drag and lift) 
and the moment of the cross sections of the girder bridges of 
standard shapes are obtained. These data refer to the span 
structures at the construction phase (without deck and parapets, 
without parapets) and operation phase, under the conditions of 
model smooth and turbulent incoming flow. The latter allows 
us to outline the boundaries for more weighted estimates of the 
aerodynamic characteristics of the girder bridges in a real wind 
current. The best agreement with the experimental data was 
obtained from the drag of the cross-section. The magnitude of 
the lifting force is more sensitive to the presence and extent of 
the separation regions, so its numerical determination is less 
accurate. The reproduction of the angle-of-attack effect on the 


aerodynamic moment of the cross-section is the most challeng- 


“The research is done within the frame of the independent R&D. 


Beedenue. BepudulupoBaHa MeTOAMKa 4YMCJICHHOrO MOJe- 
JIMpOBaHHA MOMepeyHoro OOTeKAaHHA MpOJICTHBIX CTPOCHHH 
MOCTOB Ha OCHOBe HeCcTalMOHapHoro pellileHua PevHobAca 
ua ypaBHenui Hasppe — Croxca (URANS, Unsteady 
Reynolds-averaged Navier — Stokes). JlanHHbii JByMepHbIit 
NOAXO UCHOUIb3yeTCA B COBPCMCHHBIX MeTOaXx MH MakeTax 
IIpHKIaHbIx WporpaMM BbIYHCIMTeIbHOU ruApoaspowHHa- 
MuKH. Ilenmm padoTbl — oTIaqKa WU 9KCIepHMeHTasIbHOe 
oOocHOBaHHe ykKa3aHHOH MeTOAUKH. Jia peann3aiuu m0- 
CTaBJICHHOM IesIM MUCMOUb30BaHa paHee pa3pa0oTaHHas aB- 
TOopaMu 6a3a JjaHHbIX 110 asyposmMHAMUYeCKMM XapakTepu- 
CTHKaM IMOMepeyHbIX Ce4YeHHHM MpOJICTHBIX CTpOeHHH Oa- 
JIOUHBIX MOCTOB THIOBBIX (OpM. 

Mamepuauoi u memodoi. IIpopeqeHo 4ucneHHOe MOeIMpO- 
BaHHe TOMepeyHOTO OOTeKaHHA MOCTOBbIX CTpOeHHH H3- 
KOTYpOyJICHTHbIMH (riaqKHMU) UW TypOyJeHTHBIMH BO3yII- 
HbIMM TOTOKaMM B J{Mala30He WpakKTH4YeCKH MHTePeCHbIxX 
yrioB ataku. B utTore UCHOub30BaIach MOJeb TypOyeHT- 
HocTu SST k — ow. MetToauka npezBapuTesbHO OTpaboTaHa 
Ha TeCTOBOH 3aya4e oOTeKaHHA OanOK MpAMOYIOJIbHOrO 
onmepeyHoro ceyeHusa. PacueTbI MpOBOAMIIUCh C MOMOLLbIO 
JIMI[EH3MOHHOLO MporpaMMHOro Komiiekca ANSYS. 
Pe3ylbmamel ucciedoeanus. TloKa3aHo, KaKUM OO0pa30M 
yro] aTaKH oOlpeyemseT CuIbI (MObeMHY!IO UH jOOOBOrO 
COMPOTHBJICHHA) HM MOMCHT TIOMepeyHbIxX ceyeHui Oaso4- 
HbIX MOCTOB THIIOBbIX (opm. IlouwyyeHHbIe pacueTHbIe 3a- 
BHUCHMOCTH OTHOCATCA K IIPOJICTHbIM CTPOCHHAM Ha CTayuAx 
MOHTaxa (6e€3 TIMTbI MepeKpbITHaA HU OTrpaxaeHui, 6e3 
OrpakeHHi) WH OKCIHyaTalwu B YCJIOBHAX MOJICJIbHBIX 
HaOerarollux MOTOKOB — ruaqkoro U TypOyNeHTHOrO. ITO 
MO3BOJIACT OUCPTHTb PpaHHIIbl WIA OoOee B3BEINeCHHBIX OIe- 
HOK adpOAMHAaMM4eCCKUX XapaKTePHCTHK OasIOUHbIX MOCTOB 
B PeaJIbBHOM BETPOBOM IIOTOKe. 

JIo60Boe compoTHBeHHe ceyeHHIO JeMOHCTpupyeT 
Hauslyulliee CorsacoBaHHe C OMbITHBIMUM JaHHbIMu. Besnuu- 
Ha MOWbeEMHOM CHIIbI OONee YYBCTBUTeIbHa K HaNM4HIO U 
IIPOTAKCHHOCTH OTPbIBHbIX 30H, HOITOMY ee pacueTHOe 
ollpeyemeHue MeHee TOYHO. HanOosee WpoOseMHbIM AA 
OOJIbINMHCTBa KOHQ@UrypallM ABIIACTCA BOCIIPOU3BeTeHHe 
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ing for the majority of configurations. 

Discussion and Conclusions. Comparison of the calculated and 
experimental data indicates the applicability of the URANS 
approach to the operational prediction of the aerodynamic 
characteristics of the single-beam span structures. In the case of 
multi-beam span structures, where the aerodynamic interfer- 
ence between separate girders plays an important role, the 
URANS approach must apparently give way to more accurate 
eddy-resolving methods. The results obtained can be used in 
the aerodynamic analysis of structures and in practice of the 
relevant design organizations in the field of transport construc- 


tion. 


Keywords: mechanics of fluid, gas and plasma; mathematical 
simulation; computational aerohydrodynamics, URANS ap- 
proach, bridge spans, aerodynamic characteristics. 
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BIIMAHUA yrla aTaKM Ha adpOMHAMMYeECKHUM MOMEHT ceue- 
HUA. 

OoécyacodeHue u 3aKii0o4“eHUuA. CpaBHeHWe pacueTHbIX U 
ONbITHBIX J{@HHbIX CBUeCTeJIbCTByYeT O MPpHMeHHMOCTH 
URANS-noyzxoja K OepaTHBHOMy peycKka3aHuto asporH- 
HaMM4eCKHX XapaKTepHCTHK OHOOAIOUHIX IPOJICTHbIX 
cTpoeHni. Ecau xe peub ueT O MHOTOOAaIOYHBIX Mpouer- 
HbIX CTPOCHHAX, CYIIECTBEHHY!O pOJIb UrpaeT aspoqnHaMn- 
yeckad HMHTepdepeHwMA MeXKTY OTACIbHbIMH OasKamMu. B 
9TOM Ciryyae BMecTO URANS-rogxoya cileqyeT IpHMeHATh 
Oolee TOUHbIC BUXpepa3pellaioMmlue MeTorbI. IlomyaeHHple 
pe3yIbTaTbI MOryT ObITb HCHOMb30BaHbI B UCCIeEOBAaHHAX 
a poMHaMUKH COOpyxKeHHH U B WpakTHKe MpOeCKTHBIX Op- 
raHn3alui B Cepe TpaHCHOPTHOYLO CTPOUTEIIBCTBA. 


Ku1roueBbie CJI0Ba: MexaHiKa %KUAKOCTH, ra3a UW TWJla3MBl; 
MaTeMaTHYecKoOe MOJeIMpOBaHve; BbIYMCIMTeIbHaA TUApo- 
aspoyqMHamuka; URANS-noyzxoy; mpomeTHbIe cTpoeHua MO- 
CTOB; aspoqMHaMH4eckHe XapakTepHCTHKH. 


Oopazeu ona YUumupoeanua: Toctees, FO. A. UncrenHoe 
MOJIeIMpoBaHve omepeyHoro oOOTeKaHHA  IIpOJICTHBIX 
cTpoeHuii OanouHbIx mMocToB/ TfO. A. __ Tocrees, 
A. J]. OOyxoscxui, C. JI. Canenko // Bectouk JOH. roc. TexH. 
yH-Tta. — 2018. — T.18, Ne4. — C. 362-378. 
https://doi.org/10.23947/1992-5980-2018-18-4-362-378 


Introduction. It is known that wind flow around engineering structures is, as a rule, instable turbulent in na- 


ture; different-scale eddy structures are observed in the flow [1, 2]. Nearby bluff bodies (to which, in particular, 


bridge spans belong), unsteady detached flow regions occur [3]. Accordingly, adequate modeling of the turbulence 


effects is now an important requirement for the simulation experiment techniques. 


Large Eddy Simulation, LES, and Detached Eddy Simulation, DES, are used to evaluate the aerodynamics 


of structures. However, the use of these methods is complicated by their high resource intensity, the reasons for 


which are as follows: 


— tridimensionality of the task; 


— strict requirements to the computational grid density in the near-wall region and in the “focus” region [4]; 


— restrictions on the time integration step; 


— relatively large time window length for gathering nonstationary statistics in steady state. 


At the same time, it is known [5] that for cylindrical prisms that are close in shape to beam bridge spans, the 


two-dimensional approach reproduces the basic flow properties (primary unstable mode in the body wake is essen- 
tially two-dimensional). Thus, a POD analysis (Proper Orthogonal Decomposition) of the flow near the prism with B 
/ H = 5 (H is depth of section) relative section depth was performed in [6]. As a result, 1t was established that the Ist 
and 2nd disturbance modes are two-dimensional (constant over the span) and correspond to the vorticity transfer 
along the surface. Three-dimensional modes change along the span at a reference length that is no less than B section 
depth. 

For the operational prediction of aerodynamic characteristics (ADC) of bridge structures and wind-tunnel 
test tracking, the authors used nonstationary 2D modeling based on URANS approach, Unsteady Reynolds-averaged 
Navier-Stokes. Its applicability to the definition of ADC of the bluff bodies (stationary and oscillating) was studied 
in a number of works by foreign authors (see, for example, [7]). 

Materials and Methods. When setting up computer-based experiments, the recommendations given in [8— 
10] were considered. The calculations were carried out in the ANSYS Fluent program. 

The technique was preliminary tested on the check problem of flow around beams with rectangular cross- 
section. As a result, k — @ shear stress transport (SST) model was chosen to describe the flow turbulence, and the 


grid parameters and the numerical algorithm were selected. 
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The extension of the rectangular computational domain is (30 ... 40) H lengthwise, and (14 ... 20) H — trans- 
versely. The front face of the streamlined body was spaced apart from the input boundary at (8 ... 12) H. 

We used low-Reynolds-number grids (y" < 4 ... 5 dimensionless distance of the first node to the wall) that 
enabled to calculate the boundary layer separation and reattachment. Considering the complexity of the streamlined 
body contours, multiblock grids were constructed. The internal female block consisted of quadrilateral elements 
whose density increased closer to the body surface. A layer with a structured orthogonal quadrilateral grid was gen- 
erated immediately at the wall. The wake region was covered with a grid of square cells sizing of no more than H/15 
... H/10. The cell size increased to H/4 ... H/3 to the outer boundaries. The cross-sectional perimeter contained about 
10° ... 10° cells depending on its shape. The total number of cells ranged from 40—50 thousand (for sections of simple 
shapes) to 250—300 thousand (for complex ones). An example of the computational grid near a beam of trapezoidal 


section with overlapping and fencing is shown in Fig. 1. 





Fig. 1. Example of computational grid (fragment) 


When solving the Navier — Stokes equations, the velocity — pressure relationship was implemented using the 
SIMPLE algorithm. The convection and viscous terms of the equations of flow and the transport of turbulent parame- 
ters were approximated by schemes of second-order accuracy. 

The numerical integration was carried out by an implicit time scheme of the second-order accuracy. At inte- 
gration step was (0.02... 0.04) H/V (V is incident flow velocity), 1.e., under the vortex shedding with dimensionless 
frequency f H/V = 0.1, it was approximately 250—300 times less than 1/f period, and this provided an acceptable reso- 
lution of the non-stationary flow parameters. The established vortex trail was usually formed by H/V moment (60 ... 
120). Thus, the total number of integration steps averaged 6000 ~ 10,000. To collect nonstationary statistics, a time 
interval of at least 5 periods was used. 

An example of a qualitative comparison of the computational and experimental flow patterns near the span 


is shown in Fig. 2. 
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b) 
Fig. 2. Example of flow pattern over bridge span: experiment (a), calculation (5) 


Research Results. Detailed information on ADC typical cross sections can be found in [11]. Figures 3—14 
present a comparison of the computational and experimental data on the coefficients of averaged aerodynamic forces 
(drag, lift) and the moment for some specific sections. 








Fig. 3. Coefficients of averaged aerodynamic forces. Here: Cyq is drag; Cy, 1s lifting force; C,, is moment; B and H are longitudi- 


nal and transverse section dimensions (excluding fencing); a is angle of attack 


The incident smooth flow is characterized by the intensity of 0.5%, the turbulent one — of 8%. The computa- 
tion data is represented by solid lines. 
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Fig. 4. Narrow single square girder (B/H 0.75 ratio): cross-section shape (a); smooth flow (5); turbulent flow (c) 
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Fig. 5. Narrow single square girder (B/H 1.24 ratio): cross-section shape (a); smooth flow (5); turbulent flow (c) 367 
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Fig. 6. Narrow single girder with overlapping (B/H 2.3 ratio): cross-section shape (a); smooth flow (6); turbulent flow (c) 
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Fig. 7. Narrow single girder with overlapping and fencing (B/H 2.3 ratio): cross-section shape (a); smooth flow (6); turbulent flow (c) 
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Fig. 8. Wide single square girder (B/H 1.85 ratio): cross-section shape (a); smooth flow (5); turbulent flow (c) 
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Fig. 9. Wide single girder with overlapping (B/H 3.9 ratio): cross-section shape (a); smooth flow (5); turbulent flow (c) 
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Fig. 10. Trapezoidal girder (B/H 3.09 ratio): cross-section shape (a); smooth flow (4); turbulent flow (c) 
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Fig. 11. Trapezoidal girder with overlapping (B/H 5.1 ratio): cross-section shape (a); smooth flow (6); turbulent flow (c) 
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Fig. 12. Double-girder structure (B/H 2.75 ratio): cross-section shape (a); smooth flow (6); turbulent flow (c) 
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Fig. 13. Double-girder overlapping structure (B/H 3.85): cross-section shape (a); smooth flow (6); turbulent flow (c) 
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c) 
Fig. 14. Multi-girder structure with overlapping (B/H 9): cross-sectional shape (a); smooth flow (6); turbulent flow (@) 
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Discussion and Conclusions. The analytical results show that, with some exceptions, with an increase in the 
relative width of the B / H cross section of a single-beam structure, the accuracy of the calculated prediction of its ADC 
rises. As a rule, the best agreement is indicated for the frontal resistance of the section. For most configurations, the 
computation data is slightly higher than the drag coefficient obtained experimentally. It should be clarified that for the 
considered bluff bodies, the major contribution to the cross-section drag is made by the form (pressure) drag, which is 
mainly determined by the difference in pressure forces on the upstream and leeward sides of the cross section. The ac- 
cepted theoretical approach coarsens the dynamics of the vortex structures in the zone behind the body, which leads to 
an underestimated pressure recovery in this area. 

The lift magnitude is more sensitive to the presence, extent and type (open/closed) of the detached flow regions. 
This applies especially to the span structure equipped with a slab; in this case, it is possible to re-attach the flow to the 
upper side of the slab with the formation of a closed separation zone (approximately at B/H=5). Therefore, 1n compari- 
son with frontal resistance, the calculated determination of lift force is less accurate, especially for superstructures with 
a floor slab. 

The reproduction of the angle-of-attack effect on the aerodynamic moment of the cross section is a challenge for 
most configurations. 

If the aerodynamic interference [12] occurs under the cross-flow around multi-girder spans between beams, the 
accuracy of the ADC prediction falls with an increase in the number of beams (relative overall section width). In this 
case, it is advisable to use more accurate DES and LES eddy-resolving methods instead of the URANS approach. 
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Introduction. Methods of energy saving in pneumatic drive are 
considered. The method of braking by creating back pressure 
in the exhaust cavity of the pneumatic actuator is of interest. 
Under braking, the compressed air energy is stored in the 
recuperative volume. It is possible to control the braking 
dynamics through setting the initial parameters of the 
recuperative volume. The work objective is to create a 
mathematical model describing the dynamic processes taking 
place in the pneumatic drive under braking by backpressure, 
with a constant mass enclosed in the cavities of the air motor, 
and considering variation of the initial parameters of the 
braking volume. 

Materials and Methods. A mathematical model is proposed 
that describes the speed change of the output link, pressures 
and temperatures in the cavities of the pneumatic drive 
depending on the initial parameters of the recuperative 
volume. The solution to the mathematical model is carried out 
by the numerical integration method. 

Research Results. The dependences of the output link velocity, 
pressures and temperatures in the pneumatic drive cavities on 
the initial parameters of the recuperative volume are obtained. 
Adequacy of the built mathematical model is confirmed by 
Fisher's criterion. 

Discussions and Conclusions. The results obtained can be 
used to solve the problems of energy saving in pneumatic 
drives under the organization of backpressure braking. The use 
of recuperative volume increases the technological flexibility 
of the drive during its readjustment and extends the 
possibilities of energy saving. 


Keywords: pneumatic drive, recuperation, energy saving, 
backpressure, braking, recuperative volume 
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Beedenue. PaccMoTpeHbl ciocoObl »3HeprocOepexeHHA B 
MHeEBMaTHYecKOM IpuBoye. Mutepec mpeyctapiaeT ciocoo 
TOPMOXeHHA CO3aHHeM TPOTHBOAAaBIeHHA B BbIXJIOMNHOU 
mouocTH WHeBMormpuBoya. IIpu sTOM 9HepruaA cxaTOrO 
BO3Hyxa HakallJIMBaeTCA B PpeKyliepaTHBHOM OObeme. 3ayaBasr 
HavaJIbHble MapaMeTpbI peKylepaTHBHOrO OObeMA, BO3MO2KHO 
yupaBlIaTb JMHaMMKOM TOpMoxeHHA. Leth padoTEl — 
co3qaHve MaTeMaTMYecKOM MOJeIM,  OlMChIBaIOllen 
JMHAMVYeCKHe IIpPOMecchl, IpOUcxoAlMe B IHeEBMaTHYeCKOM 
Ipupoye Uupw TOpMoxeHHH IIPOTHBOaBJIeHHeM, C 
TlOCTOAHHOU MaccoH, B3aKJIEKOUCHHOTO B TMOJIOCTAX 
MHeEBMOJIBHTaTesia BO3TyxXa, WU Cc yu4erOM MU3MeHeHHA 
HavaJIbHbIX MapaMeTPOB TOPMO3HOrLO OObema. 

Mamepuanbit u memoooi. IlpenioxeHa MaTemMaTu4ueckad 
MOJIeJIb, OIMCbIBaIOIad W3MeHeHWe CKOpOCTH JBMKeHHA 
BBIXOJHOrO 3BeHa, aBJICHHH UW TeMilepaTyp B IIOJIOCTAX 
MHEBMOIMpHBOa B 3aBHCHMOCTH OT HadaJIbHbIX WapaMeTpoB 


pekyniepaTuBHoro oObema. PemleHwe MaTeMaTH4eCKON 
MOJeJIM OCYINCCTBIIAJIOCb MeTOJIOM 4UCIIEHHOrO 
WHTerpupoBaHua. 

Pe3zyibmambi —ucciedoeanua. Ilomny4eHbl 3aBMCHMOCTH 


CKOpOCTH BBIXOJHOrO 3BeHa, aBJIeHHH WU TemMilepaTyp B 
TIOJIOCTAX MWHEBMONMpPHBOa OT HadasIbHbIX TapaMeTpoOB 
pekyniepaTuBHoro oObemMa. AJICKBATHOCTb MOJYIeCHHOU 
MaTeMaTHYeCKOH MOeM MNMOATBepxK_eHa TO KpHTepHto 
Outuepa. 

O6ocyorcoenua u 3akio4ueHnua. TlosyaeHHble pe3ysIbTaTbl MOryT 
ObITh TOJC3HbI WIA pelleHuaA 3aay 2HeprocObepexkeHuA B 
IIHEBMOIIPHBOax TIpu opraHyu3allHu TOPMO2XKeHHA 
IIpoTHBOaBeHHem. IIpumMeHeHHe pekynepaTHBHOrO OObeMa 
TOBbIMmMAaeT TeEXHOJOTM4eCKylO TMOKOCTb MpuBOa Tpu ero 


TepeHalayiKe MW paciiiupsxeT BO3MO2KHOCTH IHeprocbeperxkeHHA. 


Kyirouveeoie ciloea: WHeBMaTW4ecKHH WpuBOT, pekynepalua, 
3HeprocOepexeHHe, NpOTHBOaBIIeHHe, TOPMOxeHHe, 
pekyllepaTHBHBbIi OObeM. 
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Introduction. Pneumatic drives are widely used for automation and mechanization of secondary processes. 
Energy consumption of the pneumatic equipment may be more than 20% of the total consumption of the enterprise. 
Therefore, the issues of energy saving in pneumatic drives are urgent problems [1-3]. Such features as compressibility of 
the working environment and inertia of the output links impede shockless braking, complicate the control and drive 
design [1, 3, and 4]. It is possible to control the law of pneumatic actuator braking both by affecting a regulator device 
[1, 2, 5, and 6] and by selecting the most useful technique of braking [7-9]. 

These above features enable to use compressed air as a brake damper and to accumulate braking energy, which 
is effectively implemented under braking the pneumatic actuator through back pressure [7]. While changing the switch 
coordinate for braking, pressure in the brake and injection cavities, connecting additional volume to the exhaust cavity, 
it is possible not only to provide shockless braking, but also to recuperate the energy of compressed air under braking [2, 
7, and 10]. 

By setting the initial parameters of the regenerative volume, it is possible to affect the braking and energy-speed 
parameters of the pneumatic drive. However, this issue is not adequately investigated. 

Research objective is to create a mathematical model describing the dynamic processes occurring in the 
pneumatic drive under braking by backpressure, with a constant mass enclosed in the cavities of the air motor, and 
considering variation of the initial parameters of the braking volume. 

Problem Statement. It is required to describe mathematically the dependence of the dynamic characteristics of 
a pneumatic drive under backpressure braking on the initial parameters of the regenerative volume. 

Mathematical Dependences. Consider the backpressure braking by complete overlapping of the pump and 
exhaust lines. The disadvantage of this method is that when you change the working movement and the external load on 
the output link, you have to re-set the braking actuation coordinate. An alternative is the connection of a regenerative 
volume to the brake cavity of the pneumatic motor upon braking. It is possible to affect the final pressure in the brake 
volume and the braking path of the output link of a pneumatic drive by setting its initial parameters before braking [8]. 

The following assumptions were made for a mathematical description of the gas-dynamic processes occurring 
in the cavities of the pneumatic drive [10, 12]: pressure in the lines is assumed constant; thermodynamic processes 
occurring in the cavities of the pneumatic actuator are considered adiabatic; the working medium in the cavities of the 
pneumatic drive compressed under braking is taken as an ideal gas. 

The backpressure braking is based on the principle of creating a resistance force to the movement of the output 
link of an air motor. This is achieved by partial or complete overlapping of the pump and exhaust channels. Complete 
overlap is more efficient because there is no release of pneumatic air from the brake cavity. The dynamics of the 
pneumatic drive before braking is determined by a known system of equations that considers the subcritical and over- 
critical discharge regimes [10]. 

The pneumatic drive parameters are shown in Fig. 1. The working medium parameters are identified as follows: 
Vi» Dw I, are volume, pressure, and temperature of the air compressed in the channel, respectively; V,, p,, Ty, are 
volume, pressure, and temperature of the air compressed in the head end of the air motor; V,,, p,,, T,, are volume, 
pressure, and temperature of the air compressed in the rod or head end of the air engine; Vex, Dpexs Trex are volume, 
pressure, and temperature of the air compressed in the regenerative volume; Van, Dam, Tam are volume, pressure, and 
temperature of the free air. The air motor geometrics have the following identifiers: F’,, F’,, are head and rod useful areas 
of the air motor; f,,, f, are passage areas of the pump and exhaust lines respectively; x, xo, Xo,, are the coordinates of the 
current displacement, constant head and rod “passive” volumes, respectively; s is maximum driving stroke; P 1s 


workload. 
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Fig. 1. Pneumatic drive parameters under braking by backpressure in recuperative volume 


Under braking, there is a constant amount of the compressed air in the cavities of the pneumatic engine. We 
consider that heat exchange with the environment is insignificant, therefore, we regard the thermodynamic process to be 
adiabatic, with the adiabatic index k = /.4 [10]. For braking, the control valve spool will switch to the neutral position 
and overlap the pump and exhaust lines (f,,=0, f,=0). The switching time of the valve spool is not taken into account, the 
discharge and exhaust channels are overlapped simultaneously. The connection of the recuperative volume with the rod 
end of the air motor is considered instantaneous. 

Considering the data of the assumptions of the cavity, the working medium pressure in the head end of the air 
motor will be presented in the following form: 

Pn=Cnm!/Vn)" Prams (1.1) 
where P,,, Pnm are current and initial pressure at the time of switching to braking in the head end of the air motor; v,, Vim 
are current and initial specific volumes of the head end of the air motor; k is adiabatic index. 
Reduce equation (1.1) dividing it by the area of the air motor piston: 
Pn=(Con*%m)V/Xon +X)" P nm: (1.2) 
where Xo), Xm, X are coordinates: original, brake actuation, and current position of air motor piston, respectively. 

At the moment of shifting the distributor on braking, the regenerative and braking volumes are combined; in 

this case the compressed air parameters will be determined by the following system (1.3): 


» vk ngieil (1.3.1) 
(3) um wm 
: (13) 
P nex ” peK “e '2 


Where Vim, Vpex are specific volumes of the air compressed at the start of braking in the rod end and in the regenerative 
air volume; v;, v2 are “conditional” specific air volumes in the rod end and in the recuperative air volume; Dyim, Dpex, P 
are pressures of the air compressed at the start of braking in the rod end, recuperative volume, and in the “combined” 
volume, respectively. 

We express Vv, V2 specific volumes in the system of equations (1.3): 





k k 
k f 4 
wm wm l wm | ia 
wm wm 
| I 
k k 
=>V {> =v (p) (1.4.1) 
wm \ wm 1 
(1.4) k k 
Vy 
k k peK 9 
‘V =p:v_>p = p:| —— | > 
peK peK 2 peK | » a 
peK pexk 
1 
i k 
ay (, -v(p) 
pex \ pexk 2 (1.4.2) 
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where Vium, Vpex, are rod and recuperative volumes, respectively, at the time of the start of braking; V;, V2 are 
rod and recuperative volumes connected to the “combined” volume; m,,,,, Mex are air mass compressed in the rod and 
recuperative volume at the time of the start of braking; p is pressure in the “combined” volume. 

To determine pressure in the “combined” volume, we sum up the equations (1.4.1) and (1.4.2): 


k 
1 1 
ae k 
a (Pim) ere (P yer 
rr -_ 
1 2 


wm —peK 
fees SO cdo 1.6 





7 
wm peK 
where 7 is temperature of the air compressed under braking in the “combined” air volume; R is absolute gas constant. 
The final dependence of the air temperature in the brake chamber of the pneumatic actuator on the initial 
parameters of the regenerative volume will be obtained through determining the values of specific volumes in the 
equation (1.6) and reducing it to the area of the rod end of the pneumatic motor: 


wm ~” peK 
[P= p ((st+xp +h = (1.7) 
t peK (StXp x) T -p +h  -T <p 
PEK “um pex “™ ~ nex 


where /,,,, is reduced to the rod area of the air motor, recuperative volume; T, Tym, Tye, are absolute temperatures in the 


“combined”, rod and recuperative volumes, respectively. 














2 
d~x 
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Ke EX 
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1 u 
(s+xp5 -x)( Pum )k + Nnex (Ppex )* (1.8.3) (1.8) 
5 . . ie 
[s+299 + hpex +) 
k-1 

p k 
‘euloe 7” (1.8.4) 

P 

n 
ee 
eK 

— P((S+Xp, + nex — x) s )) (1.8.5) 


(S+Xq,, —x)T nex ‘Pu + hyex ‘Liu Ppex 


The air motor piston movement is described by the equation (1.8.1). The general system of equations for the 
braking process of a pneumatic device will have the form (1.8). 

Here: m is mass of the working body reduced to the air motor rod; 7), is absolute air temperature in the pump 
line; 7;,,, 18 initial value of the air temperature in the head end of the air motor at the time of shifting to braking. 
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The equations (1.8.2) and (1.8.3) describe the pressure variation, and the equations (1.8.4) and (1.8.5) describe 
the temperature change in the piston and brake chambers of the air motor, respectively. 

The simultaneous solution of the system of equations describing the braking dynamics of the pneumatic 
actuator (1.8) and its acceleration by Runge-Kutta numerical quadrature method [14] enables to study the dependence of 
the dynamic characteristics of the pneumatic actuator on the initial parameters of the regenerative volume. The bench 
assessment of the characteristic coincidence of theoretical and practical dependences showed satisfactory convergence 
[13], which allowed us to test the adequacy of the obtained mathematical model upon the Fisher’s variance ratio. For 
this purpose, a two-factor experiment was conducted. The initial pressure and the value of the regenerative volume were 
chosen as independent factors. The factors changed according to three levels of variation, which allowed the use of nine 
combinations of the factors with three-time replication of each experiment. 

The mathematical model is adequate to the practical results obtained through the experimental verification, 
since the calculated value of the Fisher’s criterion is 2.67, which is less than the tabulated one (2.7) [15]. 

The results obtained have enabled to proceed to a computational experiment, which will allow us to determine 
rational combinations of the initial parameters of the regenerative volume for the given law of pneumatic actuator 
braking and for maintaining high-speed parameters. 

Conclusions: 

1. A mathematical model that adequately describes the dynamic processes taking place in the pneumatic drive 
chambers under backpressure braking with accumulation of the compressed air energy into the recuperative volume is 
developed. 

2. Mathematical dependences of the pneumatic drive parameters on the initial parameters of the regenerative 
volume under backpressure braking are obtained. 
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'WHcTHTYT pawMOTeXHH4eCKHX CHCTeM H yupaBieHua KOxHoro dbexepanbHoro yHuBepcuteta, r. Taranpor, Poccuiickaa Dexepanna 


Introduction. Heuristic synthesis is used to improve the effi- 
ciency of reception and processing of discrete signals under 
aprior information pressure. The analysis of the decision- 
making algorithm for the linear-logical processing of discrete 
signals in case of the incomplete aprior data on their parame- 
ters is presented. The work objective is to develop and analyze 
the efficiency of the linear-logical algorithms. 

Materials and Methods. New mathematical algorithms for the 
signal reception and processing, effective under conditions of a 
priori uncertainty, are proposed. They are based on the consid- 
eration of the structure of emissions and process exceedance in 
the signal processing channels. 

Research Results. Linear-logical algorithms for processing 
discrete signals are developed. They are based on the consider- 
ation of one, two and more detailed characteristics of emissions 
or exceedance of random processes. 

Discussion and Conclusion. The results obtained can be useful 
in the synthesis of algorithms and devices for the signal recep- 
tion and processing. Algorithms and devices are implemented 
both in an analog form and in the form of algorithms for com- 
puters. The simulation programs for the signal processing un- 
der conditions of the considerable uncertainty of aprior infor- 
mation on the signals and the channels of their distribution are 
developed. 


Keywords: decision rules, heuristic synthesis, reception and 
processing of discrete signals, probability distribution densities, 
emissions and exceedance of random processes. 


For citation: V.S. Plaksienko. Linear-logical decision-making 
algorithm for signal processing. Vestnik of DSTU, 2018, vol. 
18, no. 4, pp. 385-391. https://do1.org/10.23947/1992-5980- 
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Beedenue. Jina ToBbiienud 3d@eKTHBHOCTH Tpvema uu OOpa- 
OOTKH JMCKpeTHBIX CHTHaJIOB B YCUOBHAX Weduuuta almpuop- 
HbIX CBeJ[eHHH TpMMeHAIOT 9BpHcTMYecKHM cuHTe3. IIpen- 
CTAaBJICH aHasIM3 aIrOpHTMa MpHHATHA pellleHnsA Mp JMHeMHO- 
OrMyecKoH OOpaooTKe JMCKpeTHbIX CHTHAaJIOB B CJly4dae, ecsIu 
allpHOpHble JaHHble OO UX NapaMeTpax Hero. Llenb pado- 
Tbl —- HOcTpoeHve HW aHamM3 93ddeKTHBHOCTH JIMHeMHO- 
JIOTH4YCCKHX AITOPHTMOB. 

Mamepuaabi u memoodoi. IIpennoxenbi 9)deKTUHBHEIe B YCIIO- 
BHAX alIPHOPHOH HeompesesIeHHOCTH MaTeMaTHUeCKHe aliro- 
PHTMbI TpwHemMa u OOpadoTKH cHrHasoB. OHH OCHOBaHbI Ha 
yueTe CTPyKTypbI BbIOPOCOB HU MpeBbILeCHHM MpOMeccoB B Ka- 
Haslax OOpaooTKu. 

Pe3zyiobmamol ucciedoeanus. Co3aHbl JIMHeMHO-sIOrM4eckHe 
asIFOPUTMBI OOPaoOTKH JMCKpeTHbIX CUrHasOB. OHH OCHOBAHBI 
Ha yueTe OHO, AByX U OOsee MeTabHbIX XapakKTepHCTHK 
BbIOPOCOB WIM MpeBbIINeHHM CyYanHbIx MpOWeccos. 
O6écyorcoeHue u 3akmo4enud. Tlosy4eHHble pe3ysbTaTbl MOTyT 
ObITb MOe3HbI B IpOlecce CHHTe3a AaJIFOPHTMOB HU YCTpOMCTB 
lIpHemMa HW OOpadoTKH curHasoB. AJIFOPHTMbI HU yCTpovcTBa 
peasIM3yeMbI Kak B aHaJIOrOBOM BHJe, Tak MU B BUA aroput- 
MOB JWId QBM. Co3azanbl TporpaMMbIl MOesMpoBaHHA pu 
OOpaboTKe CHTHaJIOB B YCNOBHAX 3HAYHTeIbHON allpvopHou 
HeonpeyesIeHHOCTH HHMopMallnu O CHrHalax UM O KaHaslax UX 
pacmpocTpaHeHna. 


Kuuroveeobie cioéa: peiiaromme lpaBusia, 3BPUCTHYeCKHH CHH- 
Te3, IIpHeM HW OOpadoTKa AHCKPeTHBIX CHTHaIOB, TJIOTHOCTU 
pactipeyesIeHus BepOATHOCTeH, BbIOPOChI U MpeBbIMIeHHA CIIy- 
YaWHBbIX IIpOMeccos. 


Oopazseu ona YUumupoeanua: [Inaxcuenko, B. C. JIuHentHo- 
NOrM4eCKHH aNrOpHTM IIPHHATHA pellleHuA pH OOpadoTKe 
curHaos / B. C. Hnakcuenxo // BectHuk JJoH. roc. TexH. YH- 
Ta. — 2018. — T.18, No4. —  C,385-391. 
https://doi.org/10.23947/1992-5980-2018-18-4-385-391. 


Introduction. Heuristic engineering synthesis of nonparametric decision rules is used to optimize sig- 


nal processing under conditions of considerable a priori uncertainty. This procedure is based on the analysis of 
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the emission parameters of random processes at the output of the demodulator of the receiver of discrete signals 
[1-7]. The emission theory 1s of considerable use in the engineering practice. However, the analysis of regulari- 
ties and detailed characteristics of emissions is a complex analytical task even at a constant or slowly varying 
threshold. In the course of the analytical approach, relations are obtained that lead to nonconvergent series, 
which explains the absence of physically meaningful results. 

At present and in the near future, the use of discrete multiposition signals is promising. They include discrete 
address systems, multiple telegraphy systems (frequency telegraphy, multiposition frequency telegraphy), and systems 
with D and E codes in which elementary parcels of 7, duration at one or different frequencies are transmitted seri- 
al/parallel in time [5, 8—10 ]. 

When receiving binary signals at the resolver input, two random processes occur, and the decision procedure is 
reduced to the problem of statistical hypothesis testing. It is necessary to determine which random process (of the com- 
pared ones) has more energy on the observation interval. In this case, the decision-making procedure can be reduced to 
comparing the difference signal value at the receiver output with a constant zero threshold. Hence, it is necessary to 
compare at least two random processes (at the best case, a random process from the output of the receiver of discrete 
signals at a slowly varying threshold is analyzed). 

Statistical testing of hypotheses is reduced to the analysis of the mutual exceedance of two or more processes. 
The analytical presentation of this problem is cumbersome, and it does not provide engineering solutions [1, 11—15]. 
Heuristic synthesis and computer-based statistical modeling enable to obtain significant engineering applications. 

When receiving multiposition signals, the following can be simplified: 

- their spectra S (f); 

- amplitude-frequency characteristics (AFC) of the filters of K (f) receivers (Fig. 1). 


K (f) 





fi f2 fm 


Fig. 1. Frequency response of filters 


A linear receiver of multiposition signals should have: 

- common intermediate frequency amplifier (IFA), 

- general decision making circuit (DMC). 

Sets of separation filters (see Fig. 1) and amplitude detectors (AD) are also required. The analysis of such re- 
ceivers shows [1] the following: the greater the number of m signal locations, the lower their immunity. Requirements 
for the frequency response of filters (see Fig. 1) are quite rigid. AFC should not be overlapped to ensure frequency or- 
thogonality. In this case, the noise at the filter outputs will be independent. There should be no overlapping regions be- 
tween the AFC filters tuned to f;, frequencies (see Fig. 1). The AFC form should have a flat area in the neighborhood of 
the resonant frequency, so that the signal spectra are not distorted. 

Some existing contradictions should be observed. Thus, narrowband filters limit the operating speed. In case of 
signal depression and the Doppler effect, the degradation of quality and even failure of communication may occur. 
Broadband filters lead to the interpenetration of the signals of the neighboring frequency channels, 1.e. the orthogonality 
is violated, and, accordingly, the reception quality decreases. 

The procedures at the filter outputs of the adjacent channels in the signal detection and processing systems are 
characterized by the statistical relationships that increase with extending the mutual overlapping of the signal spectra or 
AFC filters. A detailed examination of the frequency-sharing procedure shows that we can speak of pairwise dependent 
workflows in the frequency co-channels when processing multiposition signals in the case of the overlapping AFC fil- 
ters. This is of particular importance when used to optimize the reception of multiposition signals of nonlinear or linear 
logic procedures [1, 11-14]. 

Research Objective. Under conditions of considerable prior uncertainty, the signal frequency is known with 
limited accuracy, and the range of variation of the elementary symbol durations can change up to a thousand times. Ad- 
ditional requirement is real-time operation. For simplicity, we will restrict ourselves to the case of processing binary 
signals. 
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Decision algorithm synthesis. In papers [1, 11—13], some features of processing discrete signals under condi- 
tions of a considerable prior uncertainty of the information about a symbol duration are considered. The problem can be 
solved through the linear procedures by parallelism of the decision algorithm. The problem becomes more involved 
through increasing the required accuracy and expanding the variation range of the elementary symbol duration. 

Under such conditions, the use of the adaptive procedures requires significant amount of time to adapt; either it 
is almost unrealizable due to high prior uncertainty. Let us analyze the possibilities of the heuristic synthesis of the al- 
gorithm based on the statistical properties of the mutual exceedance of two or more random processes [1, 11—12]. 


Consider a broadband reception with integration or filtering after the detector, when PAN >>] (Af, is fi is the 


receiver bandwidth, 7 is the duration of the elementary symbol). In this case, the counts of the process at the resolver 
output can be considered near-normally distributed. The expectation and dispersion of this process are determined by 
the relations [4]: 


AM|x] = M[x,] — M[x2], D{x] = D[xi] + Dia]. 
The probability of erroneous reception is calculated from the formula [1—3, 12]: 
= =ll-(a) ; 


2§ coe , ees 
Here D(a) = ? | exp (—— Jak is Kramp function; a= is ratio of the constant component to 
m 0 


AMA 
VD J242AfT/R 


2 

; a , , ; 

the effective value of the variable; h* = oo is ratio of the signal element energy to the noise spectral density, where a 
Vv 





is the signal normalized amplitude, and v’ is the noise spectral density. 

When processing discrete signals under these conditions, only non-parametric decision algorithms can be used. 
The statistical characteristics of emissions of random processes are interesting themselves: ft, duration of the emissions 
(exceedance), 7 duration of the intervals between emissions, ¢,, values of emission maxima, etc. [3]. The optimization 
decision algorithms can be based on one of the detailed characteristics, for example: the crossing threshold number dur- 
ing the observation period, the duration of the threshold crossing intervals, etc. It is necessary to determine the informa- 
tive features of such detailed characteristics of the mutual exceedance of random processes. Technically realizable en- 
gineering metering data will optimize the algorithm for receiving and processing signals under the specified conditions. 

In the context of energy, the greatest accuracy under the conditions of considerable prior uncertainty will be 
ensured by the consideration of Sy, areas of mutual exceedance (within the energy limit). However, to implement it in 
real time algorithmically and a fortiori technically is difficult. 

The concept of emissions is a special case of the concept of mutual exceedance for two or more processes, in- 
cluding random ones. When processing binary signals, the problem of analyzing the emissions of the difference process 
with respect to the zero threshold is set. It is adequate to the task of analyzing the mutual exceedance of two envelope 
processes in the signal processing channels [3, 11]. 

To optimize the decision-making problem in real time, it is advisable to simultaneously consider two or more 
detailed characteristics of emissions, for example, the duration and exceedance value. The product of the duration by 
the level (amplitude) of the exceedance is an evaluation of the exceedance area (considering the shape factor of the ex- 
ceedance) [11]. 

We use two detailed characteristics. 

1. g relative level of excess: 

gq =X; (t)/ X;(t), where X; (¢) and_X; (4) are the analyzed processes. 

2. © relative duration of excess: 

© = ty / Tx, where Tt, is the correlation interval of the processes at the output of the linear part of the receiver 
determined by its passband. 

To consider these characteristics in combination is not so hard technically. 

Joint consideration of these characteristics is relatively not difficult. CopMecTHbIM y4eT yKa3aHHbIxX XapakTepu- 
CTHK TEXHH4CCKH CpaBHUTeJIbHO He THy eH 

At the same time, its reliability is close to the indicators obtained allowing for Sy areas of exceedance. Consid- 
ering the meaning of the term “excess”, g should be greater than 1, therefore in a binary situation, when forming two- 
parameter distributions, g is determined by the relation [11]: 
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ae x (i) ax, 7), 


q=7 20) (1) 


——3 x, (t) 2 x, (0). 


This approach allows us to analyze the multiparameter distributions of mutual excesses of two or more random 
processes. 

To improve the quality of decision-making, it is necessary to reduce the total number of mutual exceedance of 
processes. This follows from the analysis results of the two-dimensional laws of the mutual exceedance distribution of 
the signal-noise mixture envelope and noise envelope [11]. It is necessary to transform the processes in such a way that, 
without disturbing the likelihood ratio, to obtain the two-dimensional distribution forms, easily distinguished by the 
resolver [12]. 

Both problems are solved using a modified combined addition algorithm [1, 11, 14]: 


Xin (O=[X1() — K Xon OE) WA — K -X2 (OI, 
(2) 
Xon (4) = [X2(t) — K Xin] An — KX (OI, 
where K is a coefficient taking values from 0 to 1; 1 [Z(a)] is a single step function, with 1 [Z(t)] = 1 for Z() > 0 and 1 
[Z (£)] = 0 for Z (t) <0. 


The spectra of the X(t) and_X x (¢) processes formed after processing by the algorithm (2) are extended. They 
should be limited to the width of the spectra of the initial information processes X(t) and_X>(¢). 

The value of K determines the implemented modifications of the algorithm and devices of the combined addi- 
tion. We are talking about the algorithm of mutual transformation, the method of combined addition, the cross-blocking 
system [1, 11]. 

When K = 1, only the fact of the process exceedance is taken into account of all the detailed characteristics of 
the exceedance. The auto-selection algorithm is implemented when the diversity technique is used [4], and the mutual 
conversion — under decision making [11]. 

When K = 0.414, both the fact that one process is exceeded by another, and the level of g excess are considered 
[4]. The combined addition is implemented in the diversity reception. When | > K > 0, not only the fact of excess, but 
also g level of excess is taken into account. The modified method of combined addition is implemented under making 
decision. 

Linear-logical procedures implemented as a result of heuristic synthesis enable statistically to transform the 
original random processes. When K = | in (2), auto-selection of processes is realized in time. Temporally continuous 


processes (Fig. 2, a) turn into sequences of pulses with X("(t) and X$”(¢) random amplitudes (Fig. 2, b and 2, c). 
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1] 
Ka | | | 
| | | 

| | 





d) 
Fig. 2. Combined addition at K=1 


In Fig. 2, d, the difference processes are presented: Z =X, — X> without processing and Z“ = x? — x, after 
processing according to rule (2). 
Density of probability distribution of X(?(¢) and X$”(¢) processes: 


W) (x) = A.B(x) + W, (XE) (3) 
W;? (x) = A,8(x) + W,(a)F, (x). (4) 
A, and A, coefficients are determined from the normalization condition: 
A =1- i; W,(x)F, (x)dx, (5) 
A, =1-|W,(x)R(x)de, (6) 


where F/(x) and F,(x) are cumulative distribution functions. 


From the consideration of Fig. 2, a - 2, c, it follows that X(?(t) and X$?(t) processes are equal to zero for 


some time. Therefore, their W(x) and W(x) probability distribution densities will contain O(x) delta functions 


(Fig. 3). 
We 


Asc x) 





Fig. 3. Density of probability distribution of processes 


Fig. 4 shows the probability distribution density of difference processes W(z) without processing according to 


rule (2) (Fig. 4, a) and W,?(Z) after processing (Fig. 4, b). Here, A is the erasing zone of the resolver. 
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W Wt) 


a) 6) 
Fig. 4. Density of probability distribution of Z(), Z°? (1) processes 


In accordance with Fig. 2, b and 2, c, the processes after treatment according to rule (2) are impulsive. Their 
spectra will be broader than the spectra of the source processes. Thus, processing according to rule (2) does not lead to 
an increase in the quality of reception, but only improves the threshold properties [11—12, 16]. 

Their spectra will be broader than those of the source processes. Thus, processing according to rule (2) will not 
lead to an increase in the reception quality, but it will only improve the threshold properties [11—12, 16]. 

A full treatment procedure includes two operations: 

- linear-logical operation according to rule (2) (it is pertinently non-linear, since there appear components in 
the spectrum that were not at the input); 

- filtering. 

Filtering considers the second parameter of © exceedance — relative duration. The expanded spectra of pro- 
cesses remain within the boundaries specified under the formation of X(t) and_X,(t) source processes (see Fig. 2, a). 

Processing according to rule (2) can repeat depending on the variation range of the duration of the expected 
signals and reception conditions [12]. In this case, the values of K parameter are different, they are always less than 1, 
and increase in the subsequent processing cycles. 

Problems of the automatic frequency control are solved more efficiently when using devices that implement 
linear-logical processing [14, 17-18]. Locking and retention bands are extended by low-end techniques. Certain param- 
eter stability of the automatic frequency control under the additive interference is provided. 

Research Results. The results of studying the algorithm (2) in case when signals are dependent are beyond the 
scope of this paper. However, the data already obtained have shown the algorithm efficiency up to the values of the 
cross-correlation coefficients of p = 0.5—0.6 processes. This is in good agreement with the results of solving the tasks of 
the diverted reception [2-4]. Algorithms of the modified combined addition of signals are effective both under the di- 
versity technique, and in decision-making problems. Yet, the consideration of their implementation features is beyond 
the scope of the paper. 

Thus, the optimization of the decision-making procedure for processing discrete signals under conditions of 
considerable prior uncertainty can be performed on the basis of non-parametric algorithms with the heuristic considera- 
tion of the detailed characteristics of mutual exceedance of random processes. More detailed characteristics complicate 
the algorithm. Still, it is this algorithm that provides greater invariance in terms of noise immunity under varying the 
duration of the expected signals. The synthesized algorithms can be implemented in the form of additional processing 
procedures on a computer, and technically — in the form of analog devices [14, 19]. 
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Beedenue. B ctaTbe paccMaTpuBarOTCA BOIPOCbl HOPMUPOBaHHA 
WoKa3aTeyIeH HajlexKHOCTH Ha paHHHXx cTaMAX IpOeKTHpOBaHHA 
MalluH. Takow moaxoy, oOecneuuBaeT coOsoqeHve 3aqaHHOro 
YpOBHA HajlexKHOCTH pH pa3padoTKe AeTaiIei, KPHTHYHbIX C 
TOUKH 3peHHA CTaOMIbHOM 9KCHTyaTalMH. Llenb padoTbr — uc- 
CI€¢OBaTb TpoOUeMbl IpoOeKTHpOBaHHuA B ykKa3aHHOW cdepe. 
AHasIv3 Me€TOJOB HOPMHPOBaHHA HajlexKHOCTH MO3BOJIAeT YTBEep- 
*KaTb, TO UX HEJOCTATOYHO JIA WpoeKTHpoBaHUuA. Kpome Toro, 
OTMeYeHBI MPOTHBOPe4uHA, CBA3AHHBIe C MCIOJIb30BAHHEM B IIpo- 
eCKTHPOBAaHHM IKCIHOHEHIIMAJIBHOTO 3aKOHA Hajlex*KHOCTH. 
Mamepuanbi u memoovi. Ha craquu TexHuuecKkoro 3aaHHa 
OlpeyesIeHbI YMCIICHHbIe 3HAaYeCHHA HapaOOTKUH UM Oe30TKa3HOCTH 
MallIMHbI. 3aTeM 3HAYCHHA BEPOATHOCTeH Oe30TKAa3HOH padoTHI 
WIM OTKA30B CHCTeMBbI MW ee 3JIEMeCHTOB MpeCTaBJICHbI CTerIeH- 
HBIMM BbIPAKeHHAMU. B HHX CTeIeHH ABIIAIOTCA TapamMeTpamMu 
pawKUpOBaHHA MOKa3aTesieH Oe30TKa3HOCTH 3JIEMCHTOB CHCTe- 
MBI. 

Pe3yibmamol ucciedoeaHuaA. CHHTe3HpoBaH asIbTepHATHBHBIM 
TOXOJ, K HOPMHPOBaHHIO, KOTOPHIN MO3BOJIAeET BBIMOJIHUTH MOJI- 
HbIM CTPYKTYPHbI aHasu3 WpoeKta. TakuM 00pa30M, MOryT OBIT 
OI[eHeHbI WOKa3aTeIM HaylexKHOCTH BCeM CHCTeMBbI UJI JeTasen, 
Oe30TKA3HOCTb KOTOPBIX OlpeesaeT Haex*KHOCTb MallliuH. Jleta- 
IM WU [pyrve 3IEMeHTHI, OTKA3bI KOTOPBIX He BHe3alIHbl, paccMat- 
puBaroTca Oe3 MpHMeHeHHA 9KCIOHeHIMabHoro 3akoHa. IIpu 
39TOM COXpaHAeTCA CBOMCTBeEHHad eEMy IIpOCTOTa MaTeMaTHYeCKHX 
olepanui. 

O6cyorcoeHue u 3akio4eHuA. JIA TpOeKTUPOBaHHA C 3a]aHHbIM 
YPOBHeM HajlexKHOCTH HeEOCTAaTOUHO YAMCIOBOTO 3HaYeHHA Bepo- 
ATHOCTH Oe30TKa3HOH padoTs (BbP) MallIMHsI B IesIOM, IIpHHu- 
MaeMOro Ha CTaqMW TexXHHYecKOrO 3a7aHud. HeoOxo,uMbl Tpe- 
Oyemple BBP u HapadoTKH 3JIEMeCHTOB CHCTeMBI, KOTOpPbIe ABIIA- 
FOTCA HCTOUHHKAMU OTKA30B, Pa3BHBarlOIIMXCA MO pa3HbIM 3aKO- 
HaM. Ilosy4eHHbIe pe3yIbTaTbI MOTyT ObITb UCIOIb30BaHbI KaK 
IIpH MpOeKTHPOBaHHU HOBbIX MeCXaHH4eCKHX CHCTeM C 3aJjaHHbIM 
YPOBHeM HajlexKHOCTH, Tak HW pW MOJepHUu3allMn MallInH. 


Ku1r04eBble CJIOBa: HalexKHOCTb, BEPOATHOCTb, Oe€30TKA3HOCTH, 
HapaOotTka, MalliMHa, TEXHU4eCKad CHCTeMA€, 3JICMEHT. 
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Introduction. The issues on standardization of reliability indices 
in the early stages of machine design are considered. This ap- 
proach maintains the target level of reliability when developing 
parts that are critical in view of stable operation. The work ob- 
jective is to study design problems in this area. The analysis of 
the reliability regulation methods suggests that they are insuffi- 
cient to design. Besides, there are some contradictions associat- 
ed with the use of the exponential failure law in design. 
Materials and Methods. At the stage of the task order, the nu- 
merical values of the operation time and reliability of the ma- 
chine are determined. Then the values of the reliability probabil- 
ities or failures of the system and its elements are shown by 
exponential expressions. In them, degrees are the ranking pa- 
rameters of the system reliability indices. 

Research Results. The alternative approach to normalization is 
synthesized; it enables to complete a full structural analysis. 
Thus, the reliability indices of the entire system or parts, whose 
operational safety determines the machine reliability, can be 
assessed. Parts and other components, whose failures are not 
sudden, are considered without using the exponential law. This 
preserves the inherent simplicity of mathematical operations. 
Discussion and Conclusions. The numerical value of the relia- 
bility probability (RP) of the machine as a whole, taken at the 
stage of the task order, is insufficient for the reliability target 
design. The specified RP and running time of the system ele- 
ments, which are sources of failures developing according to 
different laws, are required. The results obtained can be used 
both in designing new mechanical systems with the reliability 


target, and in the modernization of machines. 


Keywords: reliability, probability, failure-free performance, 
running time, machine, technical system, element. 
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Introduction. The conditions for launching mechanical products into manufacture and the unified system for 
design documentation (ESKD) [1—2] require considering the following engineering systems (ES) reliability parameters 
under the machine designing: 

- durability - calendar time (7;)); 
- reliability - probability of failure-free operation P(A) ...); 
- lifelength t = T hours (probability of failure-free operation should be guaranteed for this period). 

The given factors characterize the machine as a whole and are coordinated at the stage of the technical specifi- 
cations (TS). The developer should affect the failsafety, which is the key feature of the reliability of nonrecoverable ES 
and forms the durability of the restorable ES during the workover intervals. 

t lifelength is the argument of P(t) reliability function. By any law of reliability, P(t) monotonously decreases 
with increasing ¢. In virtue of this, to achieve the required reliability probability (RP), the lifelength must be efficiently 
minimized and selected according to the ES behavior [3]. 

Parts and nonrepairable products, being worn out, determine the time between failures (TBF). Ideally, the TBF 
of the elements should be the same or at least be a multiple of the selected ES lifelength. In this case, T;, calendar time 
and the required T lifelength are related as follows: 

T =T,-365-K 


a I oY, (1) 
wher 7; is calendar life in years; K’.,9 is the coefficient of ES possible use p.a.; Ky, 18 the coefficient of ES use per day; 
ITB is duty factor, which is the average ratio of on-time (acceleration time and steady-state motion time) to the ES work- 
cycle time. 

High RP is possible while reducing the required lifelength. If it needs to be significantly increased, failsafety 
will have to be maintained as follows: 
- to perform maintenance more often; 
- to budget disadvantage costs for ensuring the quality of products when designing key elements [4, 5]. 

The project RP can be selected from existing industry standards, from competitive conditions, and on other 
grounds, including the conventional reliability categories adopted for the engineering products (Table 1). 
Table 1 


Reliability categories of engineering products 





When considering the reliability categories, the following groups of the ES characteristics are taken into ac- 
count [3, 6, 7]. 

1. According to the structural type of products. Technological complexes, machines, units, mechanisms and hi- 
erarchy (levels) of their assembly units and parts are considered. Under the sequential component interaction in most 
machines, the RP grows with the transition to the lower levels of the reliability structure diagram. Then the RP of the 
details under the structural analysis may approach the indices of 3-5 reliability categories. 

2. According to the ES types. The projected level of RP directly depends on the level of the manufacturer’s re- 
sponsibility before the ES consumer. Undoubtedly, the highest level of RP is planned, for example, for aircraft, chem1- 
cal machinery, medical equipment, military equipment, etc. 

3. According to types of failure effects. The projected level of RP directly depends on the potential damage 
level in case of a machine failure (economic, environmental, reputational, etc.). 

It is generally accepted [5, 6, 8-13] that the ES failures are presented as sudden. In this case, the RP is de- 
scribed by the exponential law 

P(t)=e™. (2) 
Here, A failure rate in the period of normal operation after run-in is associated with 7 mean time between system fail- 
ures 
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with failsafety parameters 


XK=(-InP(T))/T. (4) 
From here, the ES mean time between failures (MTBF) can be calculated 
T =-T/(InP(T)). (5) 


This algorithm for forming the ES reliability parameters has both protagonists and opponents [7, 10, 14]. Do 
all failures occur suddenly? How correct are the ES reliability parameters (3) - (5) arising from the exponential law? 

How correct are the ES reliability parameters (3) - (5) arising from the exponential law? These questions arise 
when determining the failure rate (4) and the corresponding ES RP (3) regardless of the estimation technique: 

- according to specified T lifelength and P(T’) RP; 

-on T MTBF; 

- according to 7 test results. 

In the literature, structural analysis of the ES reliability of the ES, which transform the input effect (IE) into 
output parameters (OP), is usually referred to as a bottom-up network analysis from the RP components to the RP sys- 
tems. Herewith, the following schematic structures are considered: sequential (Fig. 1, a), parallel (Fig. 1, b), and com- 
bined (Fig. 1, c). 


a) 
BE 
wf 








BB>Inc 1 |>|nc2 ]>...>|nci > >Inc n |>Bn 


aa FoF 


C) 
Fig. 1. Schematic structures of ES reliability under various types of component 
interaction (BB=IE; BII=OP) 


Under the sequential interaction of the ES components (see Fig. 1, a) whose failures are independent, the ex- 
ponential law (2) represents a convenient mathematical model. If P(t) system RP and components - P(t) subsystems are 


in the ratio 


PO) =| IPO. (6) 
i=l 
then the failure rates of 1 system and 4; subsystems are in the ratio 
V= >A, (7) 
i=l 


and they are ranked depending on the accepted principle of intensity distribution over the components. 

A hierarchical schematic structure developed from the nested levels of subsystems and components [3, 6, 8— 
12] enables, selectively or throughout the structure, to perform a top-down analysis [3, 9] (from ES RP to RP of compo- 
nents). The required values of the RP of components along with the required lifelength are the initial data for designing 
parts with the specified reliability level [15, 16]. 

The top-down analysis algorithm based on (2) enables to distribute the intensities and RP components over all 
nested levels of the schematic structures with equal success. This is the case for radioelectronic systems [10, 11]. But at 
the level of the mechanical system components [3], a contradiction arises: it is impossible to apply the exponential law 
if the development of degradation failure proceeds according to another law. The system failure will occur as a result of 
the sequential interaction not related to the exponential law. All the while, the development of a schematic structure for 
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the ES reliability has sense if it serves as the basis for designing parts and selecting standardized components with the 
required RP for the selected lifelength. 

When analyzing ES RP with parallel component interaction (see Fig. 1, b), the following conditions are con- 
sidered: 

- components are constantly on, 

- their failures are independent, 

- each of m components has P(t) RP, 


- each of the m components is able to accept an input effect and convert it into an output parameter of the ES 
(see Fig. 1, b). 

In this case, the ES failure will occur after the failure of the last functional component. The probability of sys- 
tem failure through the component failure probability: 


F()=[[F,(. (8) 
j=l 
From the property 
P(t)+F(t)=1 (9) 
The system survival probability through the component RP: 
Pp) =1-] [I-01 (10) 
jel 


For homogeneous components 
P(t)=1-[1-P(o]’. (11) 

If the goal is to provide the specified ES RP, and homogeneous redundant components are taken into account, 
then using the expression of the exponential law (2) gives a simple failure rate calculation of components or m number 
of components only in the ratio (11). For more complex relations (10) and combined structures (see Fig. 1, c), examples 
using the exponential law are not given. 

Simple division into components is impossible in respect to complex structures with multidirectional compo- 
nent interaction with dependent failures. Performance and reliability parameters should be determined for the system as 
a whole, and in this case, 1t may be necessary to carry out a large amount of analytics and experimental works. A graph- 
ic representation of such diagrams [7, 9] is accompanied by a reference to the complexity of the functioning model and 
the cumbersome reliability calculations (for example, a complex closed interdependent operation of IC engine units or 
of jet engine). 

Main Part. For the design ES implementation under any kind of interaction, it is necessary to produce parts 
that make up the assembly units, and then — functional modules (mechanisms, units, power boxes, control modules, 
etc.). The modern modular principle of building machines provides for modular operations: development, assembly and 
debugging, modification and modernization, repair, replacement, etc. Parts and other components are at the lower level 
of the hierarchical design system consisting of nested blocks. Using it as a reliability structure diagram, it is necessary 
to find the criterion for RP distribution over components. We should consider not P (f) function, but its numerical value 
P(T). In this case, it is possible to substitute the specified T lifelength, and to use it under the distribution of the RP nu- 
merical values over components as the replacement cost criterion for base units in the event of their failures. The re- 
placement cost may include: 

- costs of materials used for the repair, products, diagnostics and repair work; 

- failure effects estimated in money equivalent (renewal of other damaged components in case of dependent 
failures, ES idle time losses, insurance compensation for repayment, etc.). 

This study objective is to develop a top-down analysis method for reliability structure diagrams, which elimi- 
nates the above-mentioned contradictions. The ES reliability structure diagram should be a tool for the distribution of 
component RP for the specified ES RP in the initial stages of the project. Then the selected lifelength and the balanced 
component RP will be the initial data for the design with the given reliability level. For this, the structure diagram 
should meet a number of requirements. 

1. The structure diagram should be based on the construct structure. This will simplify complex functional com- 
ponent relationships. 

2. The structure diagram should contain a mathematical reliability model available for calculations in the early 

stages of the project. 
3. The reliability model of according to the structure diagram should be based not on the exponential law, but on- 
ly on the fundamental reliability properties of the technical objects and systems (6), (8) - (11). 


Machine building and machine science 


275 


http://vestnik.donstu.ru 


396 


Vestnik of Don State Technical University. 2018. Vol. 18, no. 4, pp. 392-400. ISSN 1992-5980 eISSN 1992-6006 
Becmuuk JloncKozo zocyOapcméennozo mexnuyeckozo ynueepcumema. 2018. T. 18, Ne 4. C. 392-400. ISSN 1992-5980 eISSN 1992-6006 





4. The criterion for ranking the component RP should be the cost of the component renewal in case it fails. It can 
be considered as a monetary equivalent of compensation of damage from failure. 

5. The top-down structural analysis should be applied equally efficiently for both sequentially interacting and re- 
dundant reliability diagrams, as well as for the combined ones. 

Alternative representation of numerical value of object RP. While ¢ time has P(t) value of the reliability 
function argument, it is a function decreasing by any of the known laws or obtained statically. After selecting the re- 
quired ¢ = T lifelength, RP receives P(t) = P(T) value within 0 < P(T) < 1 limits. Such a value can be represented in a 
variety of ways from which the exponential expression is selected 


P(T)=B*. (12) 
X determined by the ratio 
X = lg (P(T)) / Ig B. (13) 
Taking the base value of B = 10 degree, we obtain the expression of the numerical RP value: 
P(T) =10", (14) 
from which 
X = Ig P(T). (15) 


X exponent is called the ranking parameter of the ES RP. In the annex to the system components, X, exponents 
are also referred to as the ranking parameters of the component RP. Further in this presentation, “the simplest and most 
important case” [13] of the system reliability is considered. 

ES component RP under their sequential interaction. Consider the ES reliability structure diagram with the 
sequential component interaction (see Fig. 1, a) whose failures are independent. In this case, the system failure under 
the failure of any of n components at t = T time is expressed through the component RP according to (6): 


p(T) =P). 


In (6), product can be obtained by a variety of P(T) combinations and types of their presentation. Applying 
(14) for the numerical values of P(T) system RP and its components, P(7') = F (T)x P,(T)x...x P (7) product can be 
represented as follows: 

10° 10" x10" % x10" 
whence the connection between the system ranking parameter and the components: 
X=X,+X,+..+X,- (16) 

From the set of X, possible combinations in (16), those ones following from C, renewal costs of failed com- 
ponents are selected. At this, components with a higher renewal cost should have larger P(t) RPs. That is, at the cost of 
C,, Co, ..., C, components restoring, the set of 1/C,, 1/C,,..., 1/C, verse values should be associated with the set 
of X,, X,,.... X, components’ ranking parameters. 


The indicated ratio can be written as sums 
tL 1 
C C Cc. 
t= ; + +...4 4 ; 
-oC AG Ae 


l l l 











in which the accepted conditions for X, ranking parameter are obtained from the termwise equality of the summands: 
l 


xX Cc 


i 





X vl 
re 


The right-hand side of this equality is called the “weight coefficient of the renewal cost of sequentially interact- 





ing components”: 
a 
a= a 7 (17) 
LG 
X, ranking parameter values that meet the ranking condition: 
X, = Xa. (18) 
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The unit of measuring the renewal cost does not matter, since the cost relations are used in (17). The top-down 
analysis of the ES component RP with sequential interaction is considered in Example 1. The system presented below 
contains three components for simplicity. However, any number of components is possible for the algorithm based on 
(14-18). 

Example 1. Calculation of the ES component RP with sequential interaction (see Fig. 1, a). Initial data: 

- system RP: P(T) = 0.9; 

- number of components: n = 3; 

- renewal costs of components (in c.u.): C, = 5000, C, = 3000, C, = 2000. 

Ranking parameter for ES RP (15): X = /g P(T) = lg 0.9 = —0.04576. 

The results of the step-by-step calculation of the component RP are presented in Table 2. 

Table 2 
Calculation of ES component RP with sequential interaction 


a 0 


2000 0.0005 | 0.483871 | 0.02214 0.950297 


Cetamamd | | zara | Xa | Bx | Mella 


peal pf 0.001033 | | 0.04576 | 0 





ES component RP under their parallel interaction. Consider the following case: system crash under the 
failure of all elements with F(T) probability at t = T moment of time corresponding to the selected lifelength. Then, the 
basic property of parallel component interaction (see Fig. 1, b), according to (8), is expressed through F’(7’) compo- 


nent failure rate: 
F(T)=[|[F.@). 
j=l 


Taking by the analogy with (14) 
F(T) =10', (19) 


where Y is determined from logarithmic equation /g (F(T )) = Y /g10: 
Y=Ig F(T). (20) 
Y exponent is called the ranking parameter of the ES failure rate. When engaged the system components, Y, 
exponents are called the ranking parameters of the component failure rate. The conditions (8, 19) can be represented by 
10* =10"' x10’ x...x 10” product, whence 
Ae ae Oe tae Oe (21) 
From Y, set of possible combinations in (21), those ones following from C,, renewal costs of failed compo- 


nents are selected. At this, components with a higher renewal cost should have lower failure rate (which also means 
large RP). That is, at the cost of C,, C,,..., C,, components restoring, the set of C,,C,,...,C,, costs should be associat- 


ed with the set of Y,,Y,,...,.¥, exponents. 


The indicated ratio can be written as sums 


Y=, Y Cc C 
J=44+4+..44 wl=ear+S 


Cc, 
a a ou yas 


Here, the selected conditions for the exponents at Y, are obtained from the termwise equality of the summands: 
Y C, 


— = J . 
r >, 
The right-hand side of the equality 1s called the “weight coefficient of the renewal cost under parallel compo- 
nent interaction”: 














By as (22) 
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Y; ranking parameter values fitting the above relation: 
Y,=Yb,. (23) 
The top-down analysis of the failure rate and ES component RP with the parallel component interaction is con- 
sidered in Example 2. The system contains three components. 
Example 2. Calculation of the failure rate and ES component RP with parallel interaction (see Fig. 1, b). Initial 
data: 
- system RP: P(T) = 0.9; 
- component number: m = 3; 
- component renewal costs (in c.u.): C, = 5000, C, = 3000 u C, = 2000. 


We calculate the system reliability indices. System failure probability from (9): 
F(T) =1-P(T)=1-0.9=0.1. 
Ranking parameter values for the system failure probability from (20): 
Y =/lg F(T) =/¢ 0.1 =-1. 
The results of the step-by-step calculation are presented in Table 3. 
Table 3 


Calculation of ES component RP with parallel interaction 


22) | 103) PO) =1-E 
Se 
J 
1 


es 
5000 0316228 0.683772 
3000 0.501187 0.498813 
2000 0.630957 0 369083 

] 


—] 
Checksums and ye > »~% F(T)=|[7,@~ | Pe) =1-[ ]a-2@) 
d 
a 


Example 3. Probability calculation of the nonfailure operation of the combined ES components (see Fig. 1, c). 
Initial data: 
- number of sequentially interacting subsystems: n = 3; 


- number of parallel interacting elements of SS 3: m = 4; 

- system RP: P(T) = 0.95; 

- renewal costs of subsystems (in c.u.): C, = 5000, C, = 8000, C, = 12000; 

- renewal costs of parallel component interaction (in c.u.): C,, =3000, C,, =4000, C,, =5000, C,, 


= 6000 
We calculate the system reliability indices. Ranking parameter for ES RP from (15): 
X = lg P(T) = 1g 0.95 = —0.022276. 
The results of the step-by-step calculation are presented in Table 4. 





Table 4 
Calculation of ES component RP with sequential interaction 


er Gee Ge | em | 
CS 
Subgtem S37 [5000 0.000200 oawv7es —[-ooiosii [0975190 | 
82-8000 ——~( 0.000128 ososr22 _|-oonosi fossa | 


SS3 «| 12000 0.000083 0.204082 ~0.004546 0.989587 
Checksumsand | «dt r=Tem| 


products | «0.000408 ft sf 0.022276 | 0.95 


"<4 — subsystem selection mark for further top-down analysis. 





From Table 4, RP of the SS 3 subsystem: P3(7) = 0.989587. Failure rate of SS 3 (8): 
F(T) =1-—P,(T) =1—0.989587 = 0.010413 . The ranking parameter for the failure rate of SS 3 (22): 


Y =lg F,(T) =1g¢0.010413 = -1.982407. 
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The results of the step-by-step calculation are presented in Table 5 
Table 5 
Calculation of SS 3 component RP with parallel interaction 


2) | 103 a 
SS3 | |= 1.982407 40. 10413 0.989587 
aco — 31 [a — a 0.330401 | 0.467303 0.532697 
Component 3.2 4000 0.222222 |-0.440535 | 0.362631 0.637369 


Component 3.3 5000 0.277778 —0.550669 0.281405 0.718595 
Component 3.4 6000 0.333333 —0.660802 0.218372 0.781628 


products 


18000 1.982407 | 0.010413 0.989587 





The required lifelength of the ES and all its elements (1) should be considered, so that the component RP cal- 
culated in Tables 2—5 will be the initial data for the calculation and selection of components with the specified reliabil- 
ity level. 

It is necessary to determine the obtained resources by the corresponding reliability models with the calculated 
RPs (Tables 2, 3, 5). To do this, the dimensions and design parameters of the components (parts, standard and other 
nonrestorable items) obtained at the stage of preliminary design are considered. Then, measures can be taken to approx- 
imate the obtained and specified resources [15, 16]. 

Conclusion. The technique of top-down analysis of the machine reliability structure diagram is developed and 
tested on numerical examples. Its algorithm coincides with the sequence of design stages: from the machine RP — to the 
RP of its parts. The selection of the required lifelength and RP distribution over the parts at the initial stages of the pro- 
ject provide the design with the specified reliability level. The structure diagram can be based on the construct structure. 
This ensures a universal approach to the formation of structure reliability diagrams under various types of component 
interaction. The criterion for ranking the RP components is the component renewal in case of its failure. At any level of 
the structure diagram, it includes the monetary equivalent of the costs of materials, products, repairs, damage from fail- 
ure. They are known with sufficient accuracy at the initial stage of machine development. 

The need to use the exponential law under the normalization of reliability and the development of its structure 
diagrams is eliminated. The RP analysis relies only on the reliability fundamental properties that are common to all 
technical objects and systems. The analysis algorithm is maintained at all levels of the structure diagrams including part 
levels. 

Mathematical models of the structure diagrams are simple, available for calculations in the early project stages, 
and equally efficient for sequentially interacting, redundant and combined reliability diagrams. They provide convenient 
possibilities for algorithmization, programming, and data correction. 
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On nanoscale phenomena in the electroacoustic sputtering process" 
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K sompocy 0 Hanopa3MepHBIX ABJICHHAX B Hporecce dIeEKTpOaKyCTH4eCKOLO HallblieHusA, 
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|? JloHcKoit rocyqapcTBeHHBIii TEXHH4YeCKHM yHUBepcuTerT, r. PoctoB-Ha-Jlony, Poccniicxaa Dezepanua 


Introduction. The effect of variable parameters of the electroa- 
coustic sputtering (ELAS) process on the characteristics of the 
crystalline structure of hardening coatings is studied. The 
ELAS parameter values providing nanostructured cover coat- 
ings for machine parts and cutting tools are determined. Hard- 
ening through using such coatings allows achieving a signifi- 
cant (5-10 times) increase in the life of hardenable machine 
parts and various tools designed for mechanical processing. To 
obtain coatings with the desired properties of the surface layer, 
nanocrystalline materials should be selected. In this case, a 
certain content of the amorphous phase is permissible. 
Materials and Methods. To carry out the X-ray structural 
analysis, the X-ray diffraction Russian-made device DRON- 
3M was used. The Scherrer-Wilson method was applied to 
determine the granularity of particle blocks from the value of 
the intrinsic broadening of the diffractogram peaks. The con- 
clusions obtained in this paper are based on the method of 
separation of the affecting factor contributions into broadening 
the diffraction reflection peaks (the Warren-Averbach meth- 
od). 

Research Results. Depending on the process conditions and 
the technique for obtaining nanostructured materials, a non- 
uniquely interpretable change in the indices of the diffraction 
peaks broadening occurs, which is generally characteristic of 
nanocrystalline metals. One of the possible explanations for 
this phenomenon is the presence of a nanosized effect in the 
hardened layer. The occurrence of the nanocrystalline struc- 
ture in the sputtered layer verifies the calculated values of the 
dimensions of the coherent scattering regions (CSR). The oc- 
currence of affecting values of the misorientation angle of the 
crystal structure is verified by the CSR value for the investi- 
gated 110 and 220 reflexes, which is supported by a high per- 
centage of the amorphous phase. 


” The research is done within the frame of the independent R&D. 
™ E-mail: nb1979@mail.ru, vs_minakov@mail.ru 
*™ PaOota BbIMNOJHeHa B paMKax MHUWMaTHuBHOH HUMP. 


Beedenue. ViccneqoBano evicTBHe BapbupyeMbIX TapaMeTpoB 
Ipowecca 9IeKTpoakycTuyeckoro HanbiieHua (QJIAH) Ha 
XapaKTepHCTHKH KPHCTasJIM4eCKOH CTpyKTypbl yIIpouHsAto- 
IMX NOKpbITHH. BelaBseHbl 3HadeHHA YKa3aHHbIX WapaMerT- 
pos, oOecneuuBarollve MOJyYeHve HAHOCTPYKTYPHbIX 3all[HT- 
HbIX NOKPbITHM WeTawei MallIMH MW pexyllero MHCTpyMeHTa. 
YnpouHeHnve Tp MOMOMM MOJOOHBIX MOKPbITHM MO3BOIIAeCT 
OCTHYb 3HayMTebHOrO (B 5—10 pa3) NOBBIeHHA pecypca 
paOoTb! ympouHAeMbIX JeTaeH MW MHCTPyMeHTOB, WpeqHa3Ha- 
4CHHBIX Jit MexaHooOpaboTKH. Jia NOyAIeHHA MOKPbITHH C 
3aJJaHHbBIMM CBOMCTBaMH MOBeEPXHOCTHOrTO CJIOA CJIeAyeT BBI- 
OupaTb HaHOKpUucTaswIM4ecKHe MatTepnabl. IIpu sToM yony- 
CTHMO OlIpeyesIeHHoe COepxKaHve aMOpdHO! (a3bl. 
Mamepuaazoi u mMemoooi. Uctionb30BaH PeHTTeHOCTpyKTYpHbIi 
aHasIM3, KOTOPbIN NpOBOMIICA Ha PCeHTTeHOBCKOM udpakiu- 
OHHOM aliiilapare oOTeyecTBeHHOrO pon3BoycTBa «J[POH- 
3M». Mero Leppepa — Busicona ipuMeHed Cc IeJIbO ompe- 
J@eHUA 3CPHUCTOCTH ONOKOB YacTHI NO 3HayeHUuIO (Uu3nv4e- 
CKOrO YIIMpeHHaA MHKOB WHdpakTorpaMMpl. BpIBobI 0 pe- 
3YJIbTaTaM paOOThI OCHOBAHbI Ha MeCTOJMKe pa3jeJIeHHA BKIIa- 
OB 3HauallMx (akTOPOB B yllIMpeHve MHKOB WudpakiMoH- 
HbIX OTPaKeHHH (MeTO, YoppeHa — AsepOaxa). 

Pe3yibmamol ucciedoeaHua. B 3aBHCHMOCTH OT TeXHOJIOrM- 
4YeCKHX PeXKXMMOB HU MeTOJa MOJYYeCHHA HAaHOCTPyKTYPHBbIX 
MaTepHasIOB IIPOMCXOAUT HeEOJHO3HAYHO MHTepupeTupyeMoe 
W3MeHeHHe TOKa3zaTeslen yUIMpeHud WudpakUMOHHBIX MHMKOB, 
4TO, B OOM, XapakTepHO JIJId HAaHOKPHCTasWIM4YeCKUX Me- 
TaJUIOB. ITO ABJICHHE MO2XKHO OOBACHHTL, B YACTHOCTH, Haliu- 
4HeM HaHopa3MepHoro 3dekTa B YIIPOUHeHHOM cuioe. IIpu- 
CYTCTBHe HaHOKpHCTaIMYeCcKON CTPyKTypbI B HalriblJICHHOM 
coe HalipaMy!O NOATBepxTaeTcA BbIYHMCJICHHbIMH 3HadeHHA- 
MH pa3MepoB oOOsIacTe KOrepeHTHOrO pacceuBaHua (OKP). 
Hanmuve 3Ha4MMBIX BeIMYHH yrila pa3z30pHeHTHpOBaHHA 
CTPyKTYpbl KPHCTasWIOB NOATBepxKaeTca BesMuMHOM OKP 
JIA UCCIeOBaHHbIX pedsekcos 110 u 220, aro No_KpenssAerT- 
CA BbICOKHM IIpOIeCHTHbIM COJlep»kKaHeM aMOpMHOH (a3bl. 
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Discussion and Conclusions. The electroacoustic scattering 
method is promising for obtaining nanocrystalline structures in 
the surface and subsurface layers of the sprayed samples. The 
ELAS process variables variation leads to the parameter 
spread of the crystal lattice and coherent scattering areas. In 
this case, there is no definite trend. In the future it is expected 
to solve the given problem. First, experiments will be conduct- 
ed to determine the optimal sputtering regimes that could 


stimulate the formation of nanocrystalline structures. Second- 


Oécyacdenue u 3akiroy4eHud. MeToy 3s1eKTpOakycTHuyeckoro 
HalbUIeHHA MepCieKTMBeH JIA WOJyYYeHUA HAaHOKpUCTasIMye- 
CKHX CTPyKTyp B IIOBEPpXHOCTHOM MH IIpHMOBepXHOCTHOM CJIO- 
aX HalIbICHHbIX OOpa3zl0B. M3mMeHeHve TexHOJIOrM4ecKHx 
napamMetpoB DJIAH mpuBowuT kK pa30pocy 3HadeHH Mapa- 
MCTPOB KPHCTaJWIM4eCKON pellleTKH HU OOacTeH KOrepeHTHOrO 
paccenBaHua. B 9TOM cyIyyae He BBIABJIAeTCA ONpedesIeHHad 
TeHAecHIHMA. B WasibHevileM IpeaMouaraetca pellleHve TaHHOU 
mpoOsembl. Bo-nepBbix, OyTyT MpOBeeCHbI IKCIePUMeHTHI C 
WeHO OlpeyesweHuA ONTHMAIbHBIX pe2xKHMOB HalibluieHus, 
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CHOCOOCTBYFOIIMX # OOpa30BaHHIO HaHOKpUCTasWIM4eCcKHXx 


ly, visual observation and evaluation of the sprayed layer 
CTpykTyp. Bo-BTopbIx, MiaHupyeTca BU3yasIbHOe HadsOeHHe 
structure using electron microscopy 1s planned. M OWCHKa CIPYKTYpbI HalbWIAeMOrO COM Mp MOMOLIM 3J1eK- 


TPOHHOM MUKPOCKOMMN. 


Keywords: nanocrystalline structures, hardening coatings, 
crystal lattice, diffractometric studies, electroacoustic sputter- 


KunoueBbie = cCJIOBa: = HaHOKpHCTaswIM4ecKHe  CTPYKTYpbBI, 
YIPOUHAFOWIMe MOKPbITHA, KPMCTaM4ecKkad pellleTKa, H- 
ing. (bpakTOMeTpHyueckHe UCCIeCAOBaHHA, 93JIeKTpOaKycTHYeCcKoe 
HaliblIeHue. 
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Introduction. To obtain nanostructured hardening coatings is an urgent and long-run engineering challenge [1, 
2|. Numerous studies in the field of materials science have shown that a tangible (by several times) change in strength, 
hardness, and wear resistance of materials is possible with a decrease in crystal grains to a certain value. 

Under this investigation, it 1s expected to study the crystal microstructure of hardening coatings using the 
method of electroacoustic sputtering (ELAS) [3, 4, 5]. Such hardening involves the use of an electric spark (highly- 
concentrated power flotation) and ultrasonic longitudinally-torsional vibrations. ELAS-hardening, using a specialized 
installation, enables to increase the operational life of parts and tools by 5-10 times. 

In the paper, the ELAS operational parameters vary in order to identify their values, which contribute to ob- 
taining nanostructural protective coatings for the parts and cutting tools. To create protective coatings with the required 
properties, alloys in the amorphous state are often used. However, in certain areas, nanocrystalline materials are used 
more often than the amorphous ones. Nanocrystals relax to a far lesser extent in use of the coating, and they are highly 
competitive with amorphous materials in their properties. Obviously, to obtain coatings with the desired and stable 
properties of the surface layer, nanocrystalline materials should be selected, though with a certain content of the amor- 
phous phase. 

Materials and Methods. Resonance methods, X-ray diffraction analysis, electron microscopy, and some other 
techniques [2] are widely used to identify nanocrystalline materials in the surface layer. In the present work, the most 
accessible method based on identifying the diffraction reflection broadening under the X-ray structural analysis of sam- 
ples is used to evaluate the microstructure parameters. “DRON-3M”, the X-ray diffraction Russian-made apparatus, was 
used for X-ray diffraction analysis. The characteristic radiation of the iron anode with the release of AK, spectral line of 
the sample, as well as of the detector, was analyzed. The X-ray tube specifications were 25 kV, 5 mA. The Bragg - 
Brentano method of radiation focusing was used [6]. The selected detector parameters are as follows: motion speed is 1 
deg/min; integration time is 5s. When determining the distances between the measurement planes, the error was + 0.001 
A, which is similar to determining the position of diffraction maxima accurate to + 0.02 degrees. Samples of a cylindri- 
cal shape made of 45 steel were studied. The diameter and height of the sample was 6 mm. 

To identify samples, special marks were applied using ELAS on the front surface of each of them. The sample 
without sputtering made from 45 steel was chosen as a reference one, not the sample that is close to pure a-Fe. 

Research Results. The profiles obtained during the work are presented in Fig. 1. 
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Fig. 1. Diffraction reflections with profiles: reference with coarse grain (a); sample for sputtering (U = 16.9 V, A = 4.9 microns) (4); 
U= 13.1 V,A = 15.2 um (0); U= 12.9 V, A= 9.9 um (d); U= 12.9 V, A=5.1 um (e); U=9.1 V, A= 15.1 pm (f) 
The shape of the diffraction profiles for the reference sample is shown in Fig. 1 (a), for samples with spread 
sputtered coating on the above modes — in Fig. 1 (db) - (f). ELAS variables are electrode-anode voltage (U) and ampli- 


tude value (A) of ultrasonic frequency oscillations. 


Diffractograms of the reference and study sample show the intensity values ratio. The relationship of the crys- 
tallized and amorphous material was calculated from these data. The contribution of the texture component was not 
considered. 

Analysis of the experimental data allows us to conclude that the main phase in most of the samples under study 
is a - Fe. The Scherrer-Wilson method (the best in this case) was used to determine the grain coarseness of blocks of 
particles according to the value of the physical broadening of the diffractogram peaks [7]. 

Table 1 presents the summarized data on the research results: 

- B physical broadening, 

- d-space, 

- a dimensional parameter of the crystal lattice, 

- size of coherent scattering regions. 
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Design parameters of crystal structure and diffractogram data 
Reflex 


Table | 
110 22 
Pt | 2 | 3 | 4 5s Poe fT te | 2 | 38 | 4 |S | oe 
_ - : aie | ee 
pe ee eee 
d 
a 
| | a 


Parameter / 
il ref 
8 


0 
4 
8.87 
73 


l 4 5 l 2 5 
2.01 1.81 95 47 01 
56.916 56.86 56.926 | 56.876|56.976| 57.2 | 144.66] 144.7 | 144.78 | 144.91 | 144.53 | 144.4 
20, deg. 
,A 
2.8553 | 2.8581 | 2.8539 | 2.8567 | 2.8525 | 2.8465 | 2.8539 | 2.8561 | 2.8553 | 2.8793 | 2.8571 | 2.8571 
,A 
111 130 131 123 14 76 qi) 64 
D,nm 
Amorphous 59 55 57.5 53.5 51 59 75 79 
component 
content, % 


110 and 220 reflexes were selected as the basic ones. 

Table 1 data demonstrate significant broadening of the diffraction peaks. It is logical to assume that such a 
width in our case is, to a greater extent, due to the highly dispersed structure of the crystallites than to the dislocation 
microstresses of the surface layer. To confirm this conclusion, the contributions from these factors are estimated by the 
Warren — Averbach method [8]. The center-of-gravity shift of the peaks on the diffractograms with a simultaneous in- 
crease in a lattice parameter to 1.15—1.17% serves as an indirect proof for the above formulated hypothesis. The shift 
itself is small, but meaningful (expressed in hundredths of a degree). 

Consider the phenomena that occur when changing the values of the process conditions (electrical voltage and 
amplitude of oscillations), and their effect on the broadening value of the diffractogram peaks depending on the lattice 
dimensional parameter (Fig. 2). 
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Fig. 2. Effect of ELAS regimes variation on a lattice parameter depending on f value of physical broadening of diffractogram peaks 


As expected, a change in the deposition modes affects a lattice parameter without a definite traceable depend- 
ence. In this case, a value variations are directly proportional to the broadening value of the lines of the diffractogram 


lines. In [9], comparable proportional dependences can be observed for various oxides. 


Such ambiguously interpolated dependences are characteristic of nanostructured materials (in particular, met- 


als) and are determined by the physical method of their synthesis [1, 10] (Fig. 3). 
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Fig. 3. Effect of ELAS variables on microstructure parameters 


Perhaps, this phenomenon is explained by the occurrence of a nanoscale effect in the hardened layer. The pres- 
ence of a nanocrystalline structure in the sputtered layer is directly confirmed by the calculated values of the sizes of 
coherent scattering regions (CSR). The significant values of the disorientation angle of the crystal structure is verified 
by the CSR values for 110 studied reflex, and also, noticeably, for 220 reflex, which is supported by a high percentage 
content of the amorphous phase. Similar phenomena can be also caused by high values of disorientation angles. 

Summing up, it can be stated that all of the above phenomena (including sizes of the coherent scattering re- 
gions) directly confirm the presence of a nanoscale effect [1]. 

Discussion and Conclusions. Thus, the method of electroacoustic sputtering is promising for obtaining nano- 
crystalline structures in the surface and near-surface layers of sputtered samples. As stated in the paper, ELAS process 
parameter variation brings into existence the lattice parameter spread and coherent scattering regions without a definite 
identified tendency. In the future, this problem will be solved. This can be facilitated by the experimental determination 
of optimal sputtering regimes that ensure the formation of nanocrystalline structures. In addition, the structure of the 
sputtered layer can be observed and evaluated using electron microscopy. 
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Introduction. Rational parameters and modes of an inclined 
batch screw mixer are validated to achieve the lowest energy- 
intensive feed mixing under observance of the zootechnical 
requirements for the feed quality on uneven mixing. The 
establishment of functional dependences between parameters 
and modes enables to design power-efficient equipment for the 
on-farm feed production. 

Materials and Methods. Experimental studies of the feed 
mixing were implemented on an inclined screw batch feed 
mixer. The experimental design included variation of four 
independent factors: mixer shaft speed, filling ratio of the 
mixing chamber, mixing time, and mixing chamber angle. 
Mixing irregularity and energy intensity of the process were 
taken as optimization criteria characterizing the mixing 
efficiency. 

Research Results. The optimization criteria versus the 
variability level factor, which are two-dimensional sections of 
the second-order response surfaces, are plotted. The rational 
values at mixing irregularity of less than 5% were as follows: 
mixer shaft speed was 27.5-36.5 min-1, filling ratio of the 
mixing chamber was 0.43—0.51, mixing time was 3.0-4.2 min, 
mixing chamber angle was 22°—25°. At such parameter values, 
the mixing irregularity will be minimal, and it will be 4.10— 
4.18%, and the process intensity is from 2.08 to 2.16 kW « h/t 
Discussion and Conclusions. The dependences obtained as a 
result of the experimental studies allowed establishing the 
domain of rational design parameters and modes of an inclined 
batch screw mixer. The results obtained can be used in further 


studies under the development of initial requirements for the 


Sciences for 2013— 2020 in R&D. 


Beedenue. OOocHoBaHbl pallMOHaIbHbie WapaMeTpbl U 
P@2*KUMBI HaKJIOHHOTO WIHEKOBOTO CMeCHTeJIA MepHOAM4eCKOrO 
WeHCTBUA JIA TOCTWKeHUA HAaMMeHbIMIeH IHePrOeMKOCTU 
IIPHrOTOBJICHHA KOPMOCMeCcH TIpH ycuOBHMU coOuOTeHHA 
300TEXHH4CCKHX TpeOOBaHHH K Ka4eCTBY IIPHTOTOBJIACMBIX 
KOPMOB 10 HepaBHOMepHOCcTH CMelIMBaHHA. YcTaHOBJIeHHe 
(YHKIMOHAIbHBIX 3aBHCHMOCTeH MexKy MapamMeTpamMu wu 
P@KUMaMH TO3BOJIAeCT TpOeCKTHPOBaTb 2ZHeprosdMeKTHBHOe 
oOopy0BaHHe JJId BHYTPHXO3AMCTBCHHOrO KOMONKOpMOBOrO 
IIpOW3BOJICTBAa. 

Mamepuazol u Memoobl. DKCIIepUMeHTAJIbHbIe UCCUeAOBAaHHA 
powecca IMpHvroTOBeHHA KOPMOCMeCH ITpoOBOAMIM Ha 
HaKJIOHHOM OJIHOBaJIbHOM IIHCKOBOM CMeCHTesIe KOPMOB 
TlepHoyMyueckoro evictBua. [IaH 9KCIepHMeHTa BKIIFOUAII 
BapbHpOBaHHe 4eTbIPbMA He3aBHCHMbIMH (pakTOpaMu: 
yacTOTOH BpallleHHa Basia cCMecHuTesA, KOIdpMuIMeHTOM 
3alOJIHCHHA ©=6KaMepbl §=6CMECCHTCJIA, MpOOJDKUTeJIbHOCTBIO 
CMeIUMBaHHA WM yIJIOM HaKJIOHAa CMeCHTeJIbBHOM KaMepbi. B 
kayecTBe KpHTepveB ONTHMM3aIIMH, xapakTepH3ylIOlIux 


3) deKTHBHOCTB CMeIUMBaHHaA, ObluIM IIPHHATHI 
HepaBHOMepHOCTb CMeIMMBAHHA UM yJieJIbHad IHEPrOeCMKOCTh 
mporecca. 

Pe3yiobmamoi —=ucciedoeanus. locrpoeHbi rpaduyeckne 
3aBMCHMOCTH =6KpHTepveB- = ONITHMU3alIMH =—s«OT:~Ssé POH 
BapbipOBaHHA (bakTOpoB, TIpecTaBAroMMe cooon 


JIBYMepHblIe Ce4eCHHA MOBepXHOCTeH OTKIMKa BTOpOro 
nopayka. PatMoHasibHble 3Ha4yeHHA pH HepaBHOMepHOCTU 
CMeIlIMBaHHaA MeHee 5% COCTaBHJIM: YacTOTa BpallleHHA Balla 
— 27,5-36,5 mun |, KOIPDUUMEHT 3anOJIHeCHHA KaMepbi 
cmecuntena — 0,43—0,51, mpoaommKUTeNbHOCTb CMelIMBaHHA 
— 3,0-4,2 mMnH.; yrouI HaksIOHa CMeCHTeJIbHOM KaMepbl — 
22°—25°. [pu Takux 3HayeHuAXx WapaMeTpOB HEOJHOPOHOCTh 
CMeIIMBaHHA OyeT MHMHUMalIbHOM UW cocTaBuT 4,10—4,18%, a 
39HEproéMKOCTh porecca coctaBusa OT 2,08 qo 2,16 KBr-4/T 

O6cyocdeHue u 3akKili04eHUuA. 3aBMUCHMMOCTH, TOJy4eCHHbIe B 
pe3yibTaTe SIKCIMeCPUMeHTAaIbHBIX MCCeOBaHH, MO3BOJIMJIM 


yYCTaHOBHTb OOaCTH pallMOHaJIbHbIX KOHCTPyKTHBHbIX 
llapaMeTpoB HU pexXHMOB- HakJIOHHOrO OHOBAJIbHOTO 
INHEKOBOTO CMecuTeIA TepHoOM4eCKOrO J[eCTBHA. 


IlomyaeHHble pe3yiIbTaTbl MOTyT ObITb UCIIOJIb30BaHbI B 
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creation of new technical means with a gravitation effect of 
intensive mixing. 
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WasIbHeMIMIUX MUCCueqOBAaHHAX IIpH pa3pa0oTKe MCXOJIHBIX 
TpeOoBaHHi Ha CO3aHWe HOBbIX TeCXHHYeCCKUX CPpeJICTB C 
TpaBUTallMOHHbIM 9(—)deKTOM HHTCHCHBHOTLO CMeIIMBaHHA. 
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Introduction. An urgent task in the design of new feed mixing equipment is to ensure a lower cost of feed 
preparation [1-3]. Thus, it is necessary to reduce the energy intensity of the processes while maintaining the required 
quality of the feed. High-quality feed mixing is most effectually performed by a periodic (batch) mixer [4, 5]. Among them, 
blade mixers and screw ones [6] are widespread. This is due to the fact that the feed composition is characterized by the 
presence of dry crushed grain (more than 90%), which does not change its properties in the process of mixing with 
additives. That is, the rheological properties of the mixture can be considered unchanged throughout the experiment [7]. 

One of the advanced models of the circulation mixers is a model with an inclined bunker, inside which one or two 
screw working bodies are installed. Herewith, the reduction of the process energy intensity is achieved due to the fact that 
the forced supply of the mixed material at the upper horizontal level is replaced by its gravitation shattering from the upper 
part of the bunker [8]. Unlike horizontal mixers, the practical absence of “pressing” feed layers in the top part of the screw 
improves significantly the diffusion mixing of the product [9]. When using modern crushers [10, 11], the intensive 
movement of the material and the equalized granulometric raw material composition enables to minimize the segregation 
effects that occur under shattering [12]. 

Materials and Methods. The experimental studies were carried out on AKM-3, an inclined single-shaft screw 
batch feed mixer developed in “Donskoy” Agricultural Research Centre. AKM-3 feed mixer (Fig. 1) with the capacity of 
2.5 m’ is designed to obtain a homogeneous mixture of dry powdered components (grain, farinaceous, protein-mineral raw 
material) under the batch (periodic) operation mode. It consists of a frame, housing with a mixing chamber in which a shaft 
with a spiral tape counterflow winding 1s placed, and of blades in the top part of the bunker. It is possible to install the 
housing on the strain gauges. 





Fig. 1. General view of inclined single-screw batch mixer 
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In such a mixer, under the shaft rotation, the displacement of material particle groups from one place to another 
(convective mixing) is performed as follows: a screw transports the material to the upper part of the inclined bunker, and 
then its gravitation shattering to the lower part of the bunker takes place. Therefore, the main advantage of an inclined 
mixer, compared to vertical and horizontal ones, is the absence of dead zones between the working bodies and the walls of 
the mixing chamber [13]. 

Irregularity in the dispersion of the reference ingredient in the mixture and the process specific energy intensity were 
taken as criteria that characterize the mixing efficiency. Shredded wheat with moisture content up to 15% with dense 
loaded density of 750 kg / m® was used in the experiments as the basic ingredient (filler) of the mixture. 

As a result, a two-component mixture composed of shredded wheat - 99%, ground common salt - 1% (by weight) 
was prepared. 

Ground common salt was used as the reference ingredients to determine the mixing irregularity. Its distribution 
was determined by taking 20 samples weighing 50 g each from different points of the mixing chamber [14]. The selection 
of a specific sample was made according to GOST 13 496.0—80. The content of the dry reference component (table salt) in 
the samples was determined through the ionometric technique according to GOST 13 496.1—98. The variation coefficient of 
the actual distribution of the reference component in samples v, % (y;) was used as an index of mixing irregularity. 

Nyx (¥2) specific energy intensity of the mixing process was defined as the total energy intensity related to the 
mass of the feed mixture measured using a three-phase electric meter [15]. 

The mixer shaft speed was changed by DELTAVFd-075E frequency converter through the current frequency 
variation of the asynchronous motor. 

A three-tier second-order design was implemented during the experiment. The experiments were carried out in 
triplicate [16,17]. 

The design of the experiments involved the variation of four independent factors affecting the mixing process 
(Table 1). 

It is established that the signification of the mixing irregularity index is much higher than that of the latter factor 
which should be considered accessory. 

Table 1 


Factors and levels of variation 


Coded identification Variability range Variability interval 
Shaft speed, min | 20—40 10 


Filling factor of mixing chamber 0.4—0.6 0.1 
Mixing eyele (ime), min 
Tilt angle of mixing chamber, grade 15-35 


Research Results. After processing the experimental results, the regression equations in coded form were 





obtained, and the corresponding correlation coefficients were determined: 

y= 4.06 — 0.96x, + 0.29x, + 0.16x3 — 0.88x4 — 0.14x,x» — 0.19x,x3 +0.54x,x4 + 0.28x x3 -—0.39x 2x4 — 0.20xK3x4 
+ 0.56x"; + 0.27x79 + 0.34x73 + 0.64x°3, 

correlation coefficient R, = 0.9607; 

Vo= 2.26 + 0.99x, + 0.42x, — 0.07x3 — 0.28x4 + 0.22x x — 0.15x,x3 +0.35x,x4+ 0.13x2x3 + 0.47x2x4 — 0.08x3Xx4 
+ 0.73x°; + +0.37x72 + 0.14x73 + 0.84x°3, 

correlation coefficient R, = 0.9781. 

The model adequacy was evaluated by the Fisher criterion. The tabular value of the Fisher criterion with the 
significance level A = 0.05 equals to F = 2.8. The Fisher criterion value in the models 1s as follows: mixing irregularity - F 
= 2.51; energy intensity of the mixing process - F' = 2.65. Comparison of the calculated values with the tabulated values has 
shown that their calculated values are less than the tabulated values. Consequently, regression models adequately describe 
the process under study. The experimental values of the Cochran test do not exceed the tabulated values. Dispersions are 
homogeneous. 

Turning from the coded values of the factors (X,, X», X3, X4) to natural ones (n, k, t, g), we have obtained 
dependences of the mixing irregularity indicators (0, %) and the mixing energy intensity (NV, kW « h) on the basic factors in 
the following form: 

— mixing irregularity: 

6 = 3.76 + 11n—0.15k — 1.02t— 2.359 + 16nk — 21nt — 59ng — 0.06kt + 0.98kp + 2.35ty + 2.15n? + 4.33k° + 3.65t7 
+0.75P: 
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— mixing energy intensity: 

N=2.19 + 23n+4.13k+ 6.250 — 13nk—3.4ng —5.23kp + 0.991’ + 048K +0399’. 

The experimental data processing has resulted in plotting characteristic curves of the optimization criteria versus 
the factor variation level. They are combined two-dimensional sections of second-order response surfaces. 

Figures 3 and 4 show some two-dimensional response surfaces of the factor influence to the mixing process. In 
particular, the dependences of the effect of the shaft speed and the loading factor of the mixing chamber on the mixing 
irregularity and the process energy intensity are shown. 
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The research has resulted in determining the rational values of parameters and operation modes of the inclined 5 


single-shaft screw batch mixer with the capacity of 2.5 m° during the preparation of feed mixtures that meet the 
zootechnical requirements on mixing irregularity under the condition of minimum energy intensity of the process. The 
rational values to obtain mixing irregularity up to 5% are the following: shaft speed is 411 
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n = 27.5-36.5 min’; loading factor of the mixing chamber is k = 0.43-0.51; mixing time is t= 3—4.2 min; tilt angle 
of the mixing chamber is g = 22—25°. In this case, the energy intensity of the process was equal to 2.08 — 2.16 kWh/t. 

Besides, as a result of the experimental studies, the smallest mixing irregularity of 3.2% is set under the following 
mixing modes and parameters of the mixer: shaft speed is n = 35 min’; loading factor of the mixing chamber is k = 0.5; 
mixing time is ¢ = 3 min; tilt angle of the mixing chamber is g = 30°. In this case, the energy intensity of the process was 
2.2 kWh/ t. The smallest energy consumption of 2.1 kWh / t under mixing irregularity of 5% is obtained at the following 
values of parameters and operating modes of the mixer: shaft speed is n = 30 min’; loading factor of the mixing chamber is 
k = 0.55; mixing time is tf = 3 min; tilt angle of the mixing chamber 1s g = 25°. 

Discussion and Conclusions. The dependences resulting from the experimental studies enable to establish 
domains of rational design parameters and modes of the inclined single-shaft screw batch mixer. The results can be used for 
the development of technical means for the production of complete feed, as well as for the modernization of existing mixed 
feed industry to increase their energy efficiency. 
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Introduction. The work objective is to increase the reliability of 
the prediction methods for the lithium hydroxide behavior in the 
steam-water circuit at the thermal power plants and nuclear pow- 
er plants, and for the operational monitoring of the pH index of 
steam solutions. A method of operational control is developed on 
the basis of the conductometric measurements of the hydrogen 
index of the corrosion inhibitor vapor solutions for construction 
materials of lithtum hydroxide used at the TPP and NPS. 
Materials and Methods. A mathematical model method is used 
for the practical implementation of the high-temperature opera- 
tional control of the steam solution pH index. 

Research Results. A method for monitoring the pH of vapor 
solutions of lithium hydroxide based on the determination of 
vapor concentration through the steam condensation in the 
coolable conductivity sensor located in the vapor space of the 
steam generator 1s developed. This has significantly improved 
the accuracy of determining the lithium hydroxide concentration. 
Equations describing the change in the limiting molar equivalent 
conductance and dissociation constants of lithium hydroxide in a 
wide range of state change parameters are proposed. 

Discussion and Conclusions. The proposed on-line technique of 
testing the pH value of steam solutions, and mathematical mod- 
els for calculating the limiting molar equivalent conductance and 
dissociation constants provide an acceptable error level calcula- 
tions and the capability of measurements automation. With an 
increase in the vapor temperature up to 573.15 K, the necessity 


arises to fortify lithium hydroxide in the vapor to 10° mol / kg. 


* The research is done within the frame of independent R&D. 
™ E-mail: viadnik48@aaanet.ru 
*™ PaOota BbIMOJHeHa B paMKax MHUUMaTHBHOK HUMP. 


Beedenue. Len padoTsl — NoBbilieHve HayeéxKHOCTH MeTOJIOB 
IIpOrHoO3HpoBaHHA NOBeAeCHHA THAPOKCH Aa JIMTUA B MapoBosA- 
HOM Tpakte Ha TOC u AOC u ONnepaTHBHOTO KOHTPOJIA BOJIO- 
POHOrO NOKa3zaTesA WapOBbIx pacTBopoB. Pa3paboTaH MeTOo 
OlepaTHBHOTO KOHTpPOJIA Ha OCHOBe KOHJYKTOMeTpH4ecKHXx 
VW3MepeHHH BOJOPOAHOTO WOKa3aTe/IA WapOBbIX paCTBOPOB 
WHTHONTOpa KOPppO3HH KOHCTPYKIMMOHHBIX MaTepHasIoB Tu- 
pokcnyla JMTUA, WpuMeHsAemoro Ha TIC u ADC. 

Mamepuaaboi u memoooi. Jia WpakTH4ecKou peasu3allHu BBI- 
COKOTeMIepaTypHOro OlepaTHBHOTO KOHTPOJIA BOOPOHOrO 
WOKa3aTesIA MapOBbIX PaCTBOPOB HCIOJIb30OBaH MeTOJ, MaTeMa- 
THYeCKOFO MOJIeIMpOBaHHA. 

Pe3zyibmamo. ucciedoeanus. Pa3spad0oTaH MeTO KOHTpOIA 
BOJOPOAHOTO MOKa3aTeIA WapOBbIX paCcTBOPpOB TuypoKcuya 
JMTHA, OCHOBAHHbIN Ha OlIpesesIeHHW KOHICHTpalMu B Tape 
IyTEM KOHJeCHCAaIIMN Mapa B OxsIaxKaeMOM KOHJLYKTOMeTpH4e- 
CKOM JlaTUHKe, pa3Mel[GHHOM B TapOBOM IIpocTpaHcTBe Tapo- 
reHepatTopa, 4YTO MO3BOJMIO 3HAYMTCJIbHO MOBbICHTb TOUHOCTh 
olpeyeseHHA KOHICHTpallun rugzpoKcuya sMTuA. [Ipenox*xeHBsI 
ypaBHeHHA, ONMCbIBarOllMe H3MeHeHHe TipeesIbHOM MOJIAPHOU 
9IKBUBAJICHTHOM 3JJIEKTPOMPOBOAHOCTH WM KOHCTaHT JIMCCOLlMa- 
IMM THApOKcuya JIMTUA B WMPOKOM Malla30He U3MeHeHHA 
TlapaMeTPOB COCTOAHHA. 

O6cyacdenue u 3akiioyuenuad. IIpenqioxKeHHbIe MeTO, Olepa- 
THBHOFO KOHTPOJIA BOAOPOAHOrO WOKa3zaTesIA MapOBbIX PacTBO- 
POB MW MaTeMaTHYecKHe MOJeIM WIA pacuétTa mpesesbHOH 9K- 
BUBaJICHTHOH MOJIAPHOM 9JIEKTPOMpOBOAHOCTH WM KOHCTaHT 
WUCcOMMalnu OOecMe4UBAIOT MPHeEMJIEMY!O JIA WpakTH4ecKux 
pacuéTOB MOTpellIHOCTb HU BO3MO2XKHOCTb aBTOMATH3alIMH W3Me- 
penun. IIpu yBesmuenuu Temiepatypsbi mapa 20 573,15 K Bo3- 
HUKaeT HeEOOXOAMMOCTh B YBCJIMYCHHH KOHICHTpalluu TH ApoK- 
cua mura B Mape qo 107 Momp/kr. 
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Introduction. When organizing and maintaining water chemistry regimes (WCR) at TPPs and NPPs, one of 
the major tasks is to minimize corrosive damages to the thermal power equipment in the zone of contact with water and 
steam [1]. Addition of ZiOH in sufficient quantities into feedwater provides the creation of a hard protective lithium- 
ferritic film on the metal surface, which slows down corrosion and stress-corrosion cracking of the heat-exchange tubes 
of steam generators of NPPs with VVER-1000 [2, 3], as well as a significant increase in uptime [3]. When LiOH is add- 
ed into the feedwater of the TG — 104 drum boilers at the TPP, the formation of a protective lithium-ferritic film is also 
observed [4]. LiOH is used at the TPPs in Russia and abroad [4—6]. In advanced fourth-generation nuclear reactors — 
SCWR (supercritical water-cooled reactor), providing an increase in thermal efficiency compared to the existing ones 
(PWR) from 33 to 44%, at pressure P = 25 MPa, T temperature of water vapor is planned to be increased to 953 K [6]. 
We need data on A, limiting equivalent conductance and K, LiOH dissociation constants at high P and 7, presented in 
an easy-to-calculate form, to create on-line monitoring of the pH values of aqueous and steam LiOH solutions at TPPs 
and NPPs that determine the effectiveness of anticorrosive protection. The technique for calculating the pH of LiOH 
solutions based on the data on y electric conductivity proposed by the authors for the liquid phase [7] is unacceptable 
for the vapor phase at low p densities and C concentration due to a high error in determining y and K, [8]. 

The work objective is to increase the reliability of methods for predicting LiOH behavior in the steam/water 
circuit at TPPs and NPPs and of the operational control of the of steam solutions pH. 

To achieve this goal, equations are developed for calculating A, and K, on the saturation line and in super- 
heated steam at T up to 773.15 K; and a technique to control pH of LiOH vapor solutions is worked out. 

Materials and Methods. Compared with the traditional control associated with sampling and cooling of sam- 
ples, the operational pH control using conductometric sensors installed in the heat conductor of the steam generators 
and working at J and P close to the operating ones, is more objective and practically inertia-free [7, 8]. In pH calcula- 
tions for the liquid phase in [7], the experimental data on Kz LiOH [8], K,, ionic product [9], and ¢ dielectric capaci- 
tance for water [10] were used. C concentration of LiOH aqueous solutions was determined in terms of C LiOH=f (y 
LiOF) r=const, the dependences obtained on the basis of the experimental data on y solutions. 

Under measuring y solutions in the vapor phase (especially for small p and C), when the polar properties of the 
solvent are weakened, the dissociation degree of LiOH is small, and the resistance of the interelectrode space of the 
measure cell of the conductometric converter can become commensurate with the resistance of the electrical insulator, 
which will lead to large errors in determining y [8]. 

For the practical implementation of the high-temperature operative pH control, the necessary data on A, and 
K, can be obtained in the studied range of the state parameter variation through the mathematical simulation technique 
and be presented in the form of equations. 

Research Results. Considering the analysis of existing methods, this paper proposes a method for testing pH of 
LiOH vapor solutions based on the determination of C LiOH in steam through the vapor condensation in the conductivi- 
ty sensor with a cooled capillary [11] placed in the vapor space of the steam generator in the coolant flow. In this case, 
C LiOH in steam is equal to C LiOH in the vapor condensate and can be determined from C LiOH=f (y LiOH) r=cons; 
dependences obtained for the liquid phase [7]. Vapor condensation can significantly improve the accuracy of determin- 
ing y and C [8]. The pH calculations for LiOH vapor solutions were performed according to the procedure described in 
[12] using data from [8] on Ky. 

To assess the effectiveness of using LiOH for corrosion protection in SCWR fourth-generation atomic reactors, the 
authors [6] calculated the pH of LiOH aqueous and steam solutions at P = 25 MPa in a wide variation range of p, T, C, 
on the basis of the experimental data on K, [13-15]. Besides, an empirical equation for calculating K, molal Kon con- 
stant [16] in p variation range from 50 to 1000 kg / m’ and T from 373 to 873 K is proposed: 

log yKa = ay + ay/t ++ (a3 + ag/t) log 0d + (as + a6/t)-(logd)’, 
where 6=p/p,; T=7/T,; p., I, is the solution density and temperature at the critical point; a1, a2, 3, @4, a5, Ae are con- 
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In superheated steam, with a decrease in p and 7 K, LiOH reduces, the difference between pH H2O and pH of 
LiOH steam solutions decreases, which reduces the effectiveness of corrosion protection of the structural elements at 
high P and 7. The authors have calculated the pH values of LiOH steam solutions at P = 25 MPa; T = 673, 723 and 773 
K for C=10°, 10°, 10°, 10° and 10° mol/kg. The experimental data on K, LiOH [8], K,, 1onic product [9], and ¢ die- 
lectric capacitance of water [10] were used in the calculations. Fig. 1 shows pH dependence of aqueous and vapor LiOH 
solutions on p at P = 25 MPa for C=10 °-10’ mol/kg. 

In the area of p <330 kg/m’, the authors [13-15] did not perform experimental studies; therefore, the figure 
shows only the results of calculating pH of LiOH using the data on K, LiOH obtained by the authors’ equation [16]. 


773 re 
723 = 
873) | 673 623 3733523473373 





100 200 300 400 500 600 700 800 900 1000 
Pp, kg/m? 


Fig.1. pH dependence of LiOH aqueous and steam solutions on p density at P = 25 MPa: 1 — HO; 2,3,4,5,6 — C LiOH = 10°, 10°, 
107, 10°, 10° mol/kg; A, © is calculation of pH from data on Kg [13-15 | and [8]; 2,3,4,5,6 lines are calculation of pH according to 
Kg, data obtained from equation [16]. 


A, and Ky, values [8], which are necessary for creating equations describing the change of A, and Ky, with 
changing p and 7, are given in Tables | and 2 in comparison with the data from other authors. These values were previ- 
ously obtained by the authors on the basis of the measurements of the electrical conductivity of high-temperature LiOH 
solutions. 

To create equations describing A, behavior on the saturation line, in order to reduce the error of data approxi- 
mation, p variation range was divided into 2 parts (from 100 to 500 and from 500 to 958 kg/m’), for each of which an 
equation of the following form was proposed: 

Ay=by:p’ + by:p’ + bap + bs, 

where by, b;, b2, b3 are definable constants; p is density of the solution, kg/m? .b,, by, b2, b3 values for p = (100—500) 
and (500-958) kg/m’ are —0.000004; 0.005878; 

~2,810,820; 1942.943226 and —1.25711 10°; 0.024177411; —16.58559801; 5334,222549 respectively. For pKg=— 
logy values, the values of d,—d¢ coefficients are calculated in the equation similar to that presented in [16] and given 
above: 

pK g=d, + dy/t + (dy + dq/t)-log 198 + (ds + d6/t)-(logd)’, 

where d=p/p,; T=7/T.; p., T, are density and temperature of the solution at the critical point; d,, d2, d3, d4, ds, do are 
constants. 

When calculating d,;—d¢ constants, we used previously obtained experimental data on K, [8]. On the saturation 
line, we considered only our own experimental data. Calculations for the superheated steam at T= 673.15 and 773.15 K 
were made in two versions: in the first one, only the author's data on K, [8] were considered, in the second — all the data 


given in Table 2. 
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Table 1 
Dependence of A, LiOH limiting molal conductivity on p density 
on saturation line and in superheated steam 
4 A,:10, cm-m* / kmol 


T= 673 K 
1710[8 1675[8] 1665[8 
1675[8 1630[8] 1620[8 


[3] 

a ee 
400 [460]; 1392173] 
os [| 1452115} oe 
00 [sos ROTEL] | BIR STIS 
DTS P=238MP) | iT 
[2758] 127OTIS 
neos}P=952MPa) [Od 
OrmsP=991MPa) | SSSCSSOSC~<“C~—“‘iSSC“‘“‘CSCSSS 
Oe 1 
Tsa{I5]_@=94MPa) [| ——SS~—~SSS 
0 
a 2 

For the saturation line: d,;=119.50984; d,=—-115.33347; d3=—-174.20233; d4=167.90533; d5=-143.78368; 
d6=183.09356. 

For the first version of the calculation for the superheated steam: d,=7.6323975; d,=3.214114; d3;=—-10.9289; 
d4=2.5463453; ds=10.34132; d¢=15.190972; for the second version of the calculation: d,;=10.113404; d,=6.133638; 
d3=17.85984; d4=8.888042; ds;=8.405401; ds=14.648395. 

The percent deviations of the experimental values of A, and pK, from those calculated by the proposed equa- 
tions are presented in Table 3. 





Table 2 
Dependence of negative logarithm of Koff constant of LiOH pK, 
on p density on saturation line and in superheated steam 


oe pKq, mol/kg 
: Saturation line T= 673.15 K T=773.15 K 
100 10.04[8] 10.07[8]; 11.06 [16] 10.16[8]; 13.34[16] 
130 8.418] 8.56[8] 8.67[8] 
160 | 7 32/8] 7.698] 
200 | 6.42[8]; 6.45[16] 6.89[8]; 8.16[16] 
300 | 4.69[8]; 4.13[16] 5.20[8]; 5.64[ 16] 
379 3 A8[ 15] 
400 | 3.98[8]; 2.82[ 16]; 3.14[13] 4.31[8]; 4.01[16] 
433 3.115] 
500 : 3.38[8]; 2.81[13] 3.43[8]; 3.15[13] 
2.45[8]; 2.32[13]; 1.26[16] 2 58[8]; 2.43[13]; 1.97[16] 
650 2. 13[13] 2.25[13] 
700 2.08[8]: 2.01[8]; 1.83[13] 2.07[13] 


712 1.59/17] 
750 1.82/13] 1.89[13] 

800 1.63[ 13]; 0.343[16] 1.65[ 13]; 0.753[16] 
807 1.41[15] (P= 11.02 MPa) 


m2 | 12915] P=94MPa) | SSS 
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961 L03[15] (P=4.88MPa) | 


Table 3 
Deviation (in%) of experimental values of limiting equivalent conductance (oA) and negative logarithm of molal Koff 
constants (6pK 4) from values calculated by equations 





T=673,15 K T=773,15 K 1 T=673,15 K T=773,15 K 
Saturation line 1 calculation calculation ver- 2 calculation 2 calculation 
version sion version version 


Saturation 


| 


508] 29[8] | 


6.44[8]; 


56[16]; 
7.96[13 eae 


Se oat 

a vi 1; 8.78 a 

3. 7 5, 2A 
etet=t=tets 
16] 
J 


10.6[13 4.02[13] 





Discussion and Conclusions. As obvious from Fig. 1, pH values of LiOH obtained by the authors at 7 = 673 K 
are lower compared to those obtained through calculating [6] by a mean of 0.2; at T= 723 K, they are higher by 0.1—0.2; 
at T= 773 K, they are higher by 0.1—0.7 pH units. 

Thus, the error ratio of the data on pH of LiOH in the superheated steam obtained by the authors and presented 
in [6] goes up with increase in 7 and C. Such a mismatch may be considered satisfactory if it is remembered that with 
decreasing p to 200-100 kg/m’, the experimental data [8] error on Ky LiOH increases to 180%, and, moreover, the au- 
thors [6] did not use the experimental data on Kz equation at p <330 kg/m’. 

The deviation of the experimental values of A, and pK, from those calculated using the equations proposed in 
this paper on the saturation line does not exceed 1%, and only at p = 917 and p = 958 kg/m’, it reaches 2.27 and 1.51% 
for A, (Table 3). In the superheated steam at 7 = 673.15 and 7 = 773.15 K for the first version of the calculation using 
only the author’s experimental data on K, [8], the greatest deviation of the experimental data of pK, from the calculated 
one does not exceed 2.5% at p = 100-300 kg/m’ and goes up to 8.95% with an increase in p to 400-700 kg/m”. For the 
second version of the calculation, considering the data of all authors presented in Table 2, the greatest error of calcula- 
tion by the equation at p = 100-300 kg/m’and T = 673.15 K; T= 773.15 K reaches 15.5% for these authors [16], and at 
p = 100-300 kg/m? and T= 673.15 K; T= 773.15 K — 67.8%. 

Based on the analysis of the obtained results, it can be concluded that the use of LiOH as a corrective additive 
for pH control in order to prevent corrosion damage to structural elements of SCVR advanced atomic reactors of the 
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fourth generation creates no difficulties at 7<673.15 K and P = 25 MPa. With an increase in temperature already at T= 
773.15 K, for effective corrosion protection, it is necessary to raise C LiOH to 10° mol/kg, which puts forward the task 
to study the vapor solvability of LiOH at these state parameters. 

The proposed method of operative control of pH of LiOH vapor solutions and mathematical models for calcu- 
lating A, and K, provide the capability to determine pH using a measuring-computing complex that runs continuously 
in automatic mode; it has a short signal delay time, and it provides an acceptable for practical calculations error of 
measurement results. 

The application of such devices in modern systems of chemical-engineering monitoring at TPPs and NPPs will 
significantly improve the reliability of pH operation control and forecasting techniques of LiOH behavior of in the va- 
por-water circuit. 
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Data measurement system for non-destructive quality testing of hard alloys © 
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wy K* 
MndopMaltvouHHo-h3MepuTeJIbHbIN KOMILICKC JIN Hepaspyuraromlero KOHTPOJIA KavecTBa TBEPAbIX CILIaBOB 


VW. K. Upi6pnii', B. JI. Basmxos’, B. U1. Urnatenko’ 
1,2,3 
Introduction. It 1s known that the quality of products from 
sintered hardmetals, if the fabrication technique is not 
violated, is determined by the phase composition and an 
average grain size of the carbide phase. However, hard alloys 
have a disadvantage inherent in all products of powder 
metallurgy - the inhomogeneity of the structure and the 
corresponding variation of properties. Traditional methods of 
monitoring the structure and phase composition according to 
the results of selective destructive tests prevent from receiving 
quality data under the production conditions and do not 
guarantee the identity of the properties of all batch products 
under study. The major method of non-destructive quality 
control of hard alloys is coercimetry, but domestic 
coercimeters are currently not produced. In this regard, the 
work objective is to create a domestic data measurement 
system for non-destructive quality testing of hard alloys and an 
assessment of the reliability of the results obtained with its 
help. 

Materials and Methods. Cylindrical and spherangular rods 
with the diameter of 4 to 10 mm, made of alloys of VK8 and 
VKIOHOM grades, were used in the work. The method of 
comparison with a certified product was used to determine the 
coercitive force in the data measuring system. 

Research Results. The coercive force of two batches of 
products from hard alloys of VK1OHOM and VK8 grades was 
measured using the developed data measuring system and 
KOERZIMAT 1.097 HcJ coersimeter. 

Discussion and Conclusions. The experimental studies show 
that the values of the coercitive force of the samples obtained 


on different equipment have approximately the same level. 


" PaOota BBITOHeHa B paMKax HHHUMaTHBHON HUMP. 
“ The research is done within the frame of independent R&D. 
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JloHcKonw rocyqapcTBeHHBIM TeXHHY4eECKHM yHuBepcute;r, r. Poctos-Ha-Jlony, Poccuiickaa Denepana 


Beedenue. Vi3BecTHO, YTO Ka4eCTBO H3JeMN M3 CIeC4YCHHbIX 
TBepJIbIX CIaBOB, IPH OTCYTCTBHH HapylleHuH TexXHOJIOrMU 
W3POTOBJICHHA, OIIpeeAeTCA (a30BbIM COCTaBOM UH CpesHHM 
pa3MepoM 3epeH KapOuysHON (ba3bl. OWHaKO TBepIble CIIaBbI 
oOaaloT HEOCTaTKOM, MpHcCylMM B TOM WIM HHOM cTeneHu 


BCeM Y3,eIMAM MOpOuIKOBOM 


CTpyKTypbI 4 
BapMalwew cCBOMCTB. TpaMI[MOHHbIe MeTOAbI KOHTPOIA 
CTpyKTypbI u da30Boro 
BbIOOPOUHBIX pa3spyllalOWIMX UCIbITAHHH, He MO3BOJIAFOT 


MeTaJIIypruu — 
HeOJJHOPOHOCTBIO COOTBETCTBYFOINeH 


cocTaBa, 0 pe3yjIbTaTaM 


ollepaTuBHO  B  YCJIOBHAX  MPOH3BOJICTBAa MOy4aTb 
WHOpMalluMto O KayecTBe HM He rapaHTUpylOT UAeCHTHYHOCTH 
CBOMCTB BCeX W3euMH UccleqyemMow WapTuu. OcHOBHBIM 
MeTOJOM Hepa3pyllarollero KOHTpPOJIA KayecTBa TBepJIbIxX 
CIIaBOB ABJIACTCA U3MeEPeHHe KOSPUMTHBHON cCuJIbI. OnHaKo 
OTCYECTBCHHBIX KOSPIMTHMeTpOB B HacToslee BpeMaA He 
mpou3Boyntca. B cBa3H C 9THM I[eJIbIO HaCTOALIeH padorTsl 
ABJIACTCA CO3aHve OTeYeCTBeHHOrO UHOpMallMOHHO- 
YW3MepUTeJIBHOIO KOMIMIeKCa JIA Hepa3pyllarollero KOHTPOIA 
KauecTBa H3]{eIMM U3 CICUCHHBIX TBEPbIX CIJIaBOB WU ONeHKa 
JJOCTOBEpHOCTH MOJIYYCHHBIX C e€rO MOMOLIbIO pe3yJIbTaTOB. 
Mamepuanbi u  memodvl. B  padote UcnOmb30BaHbI 
IWIMHIpHyeckue U ciepoKOHHYecKHe CTepxKHH JMaMeTpOM 
oT 4 yo 10 MM, M3roTOBJIeHHbIe 43 CMIaBOB MapoK BK8 u 
BKIOXOM. Ia ompezeneHua KOSPIMTHBHOM CHIBI B 
WHOpMalMOHHO-H3MepHTeIbHOM KOMIIJICKCe UCMOIb30BaH 
MeTOJ CbaBHeHHA C ATTECTOBAHHBIM H3eJIMeM. 

Pe3yibmamoi ~—ucciedoeanua. IlpopeneHbl u3MepeHua 
KOSPUMTHBHOUM CHJIbl AByX WapTHv W3qeuH U3 TBepybIxX 
cilaBoB Mapkw BKIOXOM u mapxu BK8, BpmmouHeHHBIX C 
TIOMOLLbIO pa3pa0oTaHHoro YHOpMallMOHHO- 
W3MepUTeJIBHOIO KOMIMIeKca UW KOspuNTUMeTpa KOERZIMAT 
1.097 HceJ. 
O6cyncdenue u 3QK0UeHUA. IIpopeyeHupte 
9KCIePHUMeHTAaJIbHble UCCIeEOBAHHA WOKa3asIM, 4TO 3Ha4eHHA 
KOSPUMTHBHOU CHJIbI OOpa3il0OB, MOJyYeHHbIe Ha pa3HOM 


oOopyOBaHHH, HMCIOT IIPHMepHO OMHAKOBBIM YPOBeCHb. 
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The inhomogeneity of the structure characteristic of sintered 
hardmetals being a fabrication technique effect is confirmed. 
Providing that the relative error of measuring the coercitive 
force for each sample should be within + 6%, the results 
obtained can be considered satisfactory. It 1s shown that the 
measurement algorithm proposed and implemented in the data 
measurement system through comparing the characteristics of 
the tested and certified samples with the known coercitive 
force, allows the manufacturer of hard alloy products to 
expand the sample of inspected products up to 100%, and to 
significantly reduce the costs for non-destructive testing of 
products. 


Keywords: hard alloys, non-destructive quality control, 
coercitive force, 
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TloaTBep»xeHa, KaK CJICHCTBUC TeCXHOJIOPTHH NpOU3SBOJICTBa, 


XapakTepHaa A CIICYeHHBIX TBEPJIBIX — CIIABOB 
HeEOJHOPOHOCTb CTpyKTypbl. YUnTHIBAad, UTO OTHOCHTEeIIbHAaA 
TOrpelHOCTh W3MepeHHH KOSPIMTHBHOM CHIIbI WIA KaxKOrO 
oOpa3ila OJDKHa HaxOMTbCA B TMpememax + 6 %, MO%KHO 
C4MTATb TIOJYYCHHbIe pe3yIbTaTbl YOBJICTBOPUTeJIbHBIMH. 
TIpeMIOX#KeCHHBIN WU pealIv30BaHHbIM B 


ayIrOpHTM 


Iloka3aHo, 4TO 
WHDOpMalMOHHO-H3MepHTesIbHOM  KOMIIJICKCe 
W3MepeHHaA ITyTeM CpaBHeEHHA XapaKTepHCTHK UCIIbITYeMBIX U 
aTTECTOBAHHOrO OOpa3ila C H3BECTHOM KOIPIMTHBHON CuJION, 
TO3BOJIAeCT MpOW3BOUMTeIO W3JeIMM U3 TBepybIX CIWIaBOB 
PaclMpuTb BbIOOPKy KOHTPOJMpyeMbIX V3euHH BIIOTb 0 
100% u cyljecTBeHHO CHH3HTb 3aTpaTbl Ha Hepa3pyllarolinu 


KOHTpOJIb MPOYKUMH. 


KoroueBble CJI0Ba: TBepAble CiaBbl, Hepa3pyllaroluiui 
KOHTPOJIb KayeCTBa, KOSPUMTHBHaA CHJIa, HHPOpMallMOHHO- 
M3MepHTeJIBHbIM KOMIIJICKC. 


Oopaszeyu ona yYUumupoeanua: Lpi0pui, U. K. Uxdopmatnon- 


HO-H3MepHTeJIbHBIM = KOMIUIeKC JIA Hepa3pyliarollero 
KOHTpOJIA Ka4yecTBa TBepyabIx cnyaBoB / VU. K. Lp10pui, 
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Introduction. Among the materials obtained through powder metallurgical techniques, hard alloys are widely 
used. Their basic composition is finely-divided tungsten and titanium carbides on cobalt binder. The use of products 
from hard alloys is due to a number of valuable properties, such as hardness, resistance to friction wear, wearability. At 
the same time, the breakdown in the production of powders and sintering features causes the occurrence of structural 
inhomogeneity within one batch of finished products, and, as a result, heterogeneity of the long-range operation factors. 
Consequently, sampling check of several hard-metal products from a batch using, for example, metallographic analysis 
or mechanical testing does not guarantee the identity of the structure and properties of the entire batch. 

It is known that the quality of products from sintered hard alloys is mainly determined by the phase 
composition, 1.e. the ratio of the amount of the ferromagnetic cobalt phase and paramagnetic carbides, and the average 
grain size of the carbide phase [1-3]. At the fixed content of a cobalt phase for each grade of carbide material, the 
finished product strength is almost completely determined by the carbide graininess [4]. 

An alternative to destructive testing is non-destructive quality control of hard alloys in terms of coercive force, 
which at a first approximation inversely depends on the carbide grain fineness. The coercive force value, along with the 
strength characteristics of hard alloys, is introduced into the current standards [5, 6]. However, the absence of 
production of domestic coercimeters forces makes manufacturers of hard alloy products use foreign analogues for non- 
destructive testing. 

Accordingly, the work objective is to create a data measuring system for non-destructive quality control of 
products made of sintered hard alloys, and the validation of the control results obtained with it through comparing with 
the results obtained using KOERZIMAT 1.097 HcJ coercimeter from “Institut Dr. Foerster GmbH & Co. KG” [7] 
which is included in the Russian National Register of Measuring Equipments [8]. The results of the coercive force 
measurement on KOERZIMAT 1.097 HcJ coercimeter are provided by Serpukhov Tool Plant “T'VINTOS”. 

Materials and Methods. Test items were cylindrical and spherangular rods with the diameter of 4 to 10 mm 
made of alloys of VK8 and VKIOHOM grades. 

Figure 1 shows the schematic structure of the developed data measurement system [9, 10] which consists of six 
functional blocks. 


Tsybriy I. K. and the others. Data measurement system for non-destructive quality testing of hard alloys 
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I . . . 
—~ pe — faceplate indication 





Fig. 1. Schematic structure of data measurement system 


The stabilivolt is assembled on a step-up transformer and is designed to increase the supply voltage up to 300V 
and to feed it to the magnetization coil of the measuring unit. The transformer is turned on by a solenoid starter with a 
quench circuit, which ensures cushion make-break of the measuring unit. The magnetic starter is controlled by a 
microcontroller through a triac connection circuit. 

The gage block is a solenoid with a magnetization coil and a measuring circuit fixed inside it (Fig. 2). The 
magnetization coil parameters enable to obtain the magnetic intensity inside the solenoid up to 350 kA/m. The magnetic 
strength generated in the magnetization unit solenoid enables to measure the coercitive force close to the saturation field 
of the material. 

The scaler is designed to adjust the measuring signal level for the capabilities of the microcontroller. The 
power supplier and relay automation unit is designed to input various types and levels of voltages to the components of 
scalers, mathematical processing units, and control and indication elements. 


Magnetization coil Control product 











rs 


TR OT LO 


QO 


Test coil | Balance coil 


Fig. 2. Measuring unit 


The mathematical processing device is designed to calculate and store the signal level corresponding to the 
coercitive force of the certified sample and to determine the signal level corresponding to the coercitive force of the 
tested product through the comparison technique. 

The data measurement system algorithm is as follows: 

- certified sample is placed in the measuring unit, the coercive force value is displayed on the indicator 
faceplate; 

- U, voltage corresponding to the known coercitive force value is measured and stored; 

- certified sample is removed from the measuring unit; 

- test sample is placed in the measuring unit; 

- U> voltage corresponding to the coercitive force value of the test sample is measured and stored; 

- coercitive force value of the test sample obtained through comparing U, and U> voltages is displayed on the 
indicator faceplate; 

- test sample is removed from the measuring unit. 

If it is necessary to measure a batch of products of the same grade and of one size with the certified sample, the 
measurement process is repeated starting with the 4" point. 

Research Results. The comparative analysis results obtained through measuring the coercitive force values of 
24 product samples from hard alloys of VKIOHOM grade using the developed data-measuring system (row 1) and 
KOERZIMAT 1.097 HcJ coercimeter (row 2) are presented in Fig. 3. Similar results for 46 product samples from hard 
alloys of VK8 grade are presented in Fig. 4. 
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Fig. 4. Coercitive force of samples of VK8 grade 


Discussion. The experimental studies have shown that the coercitive force values of the samples obtained 
using the KOERZIMAT 1.097 HcJ tool and the data-measuring system have approximately the same level. 

Taking into account that the relative error in measuring the coercitive force for the KOERZIMAT 1.097 HcJ 
coercimeter is within + 6%, the comparative tests results can be considered satisfactory. 

The variation of coercitive force values of samples within one grade up to 10-12% can be caused by the 
following reasons: 

- use of powder mixtures from different manufacturers under sintering; 

- ignorance of the required process conditions; 

- effect of sample configuration on the instrumental data, since lower of the coercitive force values were observed 
on spherangular samples as a result of a larger demagnetization factor characteristic of this type of samples. 

Conclusion. Thus, comparative tests of the data-measuring system have shown the possibility of its use for 
quality control of products from sintered hard alloys. The proposed integrated measurement algorithm is implemented 
through the comparison of characteristics of the test samples and the certified sample with the known coercitive force. 
This allows the manufacturer of hard alloy products to expand the sample of controlled products up to 100%, to reduce 
time on control operations, and to significantly reduce the costs on nondestructive product inspection. 
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Introduction. The theoretical description of the ion transport in 
membrane systems in the galvanostatic mode is presented. A 
desalting channel of the electrodialysis apparatus is considered 
as a membrane system. The work objectives are the develop- 
ment and verification of a two-dimensional mathematical 
model of the stationary transport of salt ions in the desalting 
channel of the electrodialysis apparatus for the galvanostatic 
mode. 

Materials and Methods. A new model of ion transfer is pro- 
posed. It is based on the Nernst — Planck — Poisson equations 
for the electric potential and on the equation for the electric 
current stream function. A numerical solution to the boundary 
value model problem by the finite element method is obtained 
using the Comsol Multiphysics software package. 

Research Results. The developed mathematical model enables 
to describe the stationary transfer of binary salt ions in the 
desalting channel of the electrodialysis apparatus. Herewith, 
the violation of the solution electroneutrality and the formation 
of the dilated domain of space charge at overlimiting currents 
in the galvanostatic mode are considered. A good agreement 
between the physicochemical characteristics of the transfer 
calculated by the models for the galvanostatic and potenti- 
ostatic modes implies adequacy of the constructed model. 
Discussion and Conclusions. The developed model can inter- 
pret the experimental study results of ion transfer in membrane 
systems if this process takes place in the galvanostatic mode. 
Some electrokinetic processes are associated with the appear- 
ance of a dilated domain of space charge at overlimiting cur- 
rents. When describing the formation of this domain, it is pos- 
sible to find out how the processes dependent on it affect the 


ion transfer in the galvanostatic mode. 


Beedenue. CraTbA TOcBALIeHa TeOpeTHYeCKOMYy OIMCaHHlo 
Impoljecca MepeHoca HOHOB B MeMOpaHHbIX CHCTeMAaX B TaJib- 
BaHOCTaTHYeCKOM pexume. B kayecTBe MeMOpaHHOM CHCTeMBI 
paccMaTpuBaeTca KaHall OOeCCOIMBAHHA 9BJIEKTPOAMAaIM3HOrO 
annapata. Llemu paOoTsi: co3qaHve U Bepudukallua JByMep- 
HOM MaTeMaTHM4eCKOH MOJeIM cTalMoHapHoro mepeHoca 
WOHOB COJIM B KaHasie OOeCCOMBAaHHA IICKTPOAMAIM3HOrO 
allapata JIA rasIbBBaHOCTaTHYeCKOrO perxkKuMa. 

Mamepuaabi u memoooi. TpeqnoxeHa HOBAA MOJIeJIb TepeHo- 
Ca WOHOB. Ona OCHOBaHa Ha CHCTeMe ypaBHeHun HepHcta — 
IInanHxa — Ilyaccona fia seKTpuueckoro MOTeHIIMasla UW Ha 
ypaBHeHun Jia (YHKUMM 3sIeKTpuyeckoro ToKa. IlomyaeHo 
YMCIICHHOe pellieHve KpaeBOM 3aa4yH MOC MeTOOM KO- 
HC4HBIX 3JIEMCHTOB C MOMOINbIO TpoOrpaMMHoro fakeTa 
Comsol Multiphysics. 

Pe3zyiobmamoi ucciedoeanua. Pa3spadoTaHHad MaTeMaTH4yecKkad 
MOJe@JIb MO3BOJIACT ONMCaTb CTal[MOHapHbIN MepeHOC MOHOB 
OUHapHOH COJIM B KaHasie OOeCCOIMBAHHA IICKTPOTMaIM3HOTO 
anmapata. IIpu 3TOM y4UTHIBaIOTCA HapylleHve 3JIEKTpOHen- 
TpasIbHOCTH pacTBOpa U (POpMNpOBaHHe paciiiMpeHHou oOsa- 
CTH IIpOCTpaHCTBeHHOLO 3apsjja pW CBepxMpeeIbHbIX TOKAX 
B rajIbBaHOCTaTH4eCKOM pexuMe. OO aleKBaTHOCTH TOCTpo- 
CHHOM MOJIeJIM CBUJETeEIbCTBYeT XOpolee CoBMayeHve Pu3u- 
KO-XHMMYeCCKHX XapakTepHCTHK TepeHoca, paCCYMTaHHBbIXx 10 
MOJICAM JIJIA TasIbBaHOCTaTHYeCKOLO HW MOTCHIMOCTAaTHYeCKO- 
TO pexKHMOB. 

O@cyacdenue u 3aKroueHUA. Pa3padboTaHHad MOJICJIb M03B0- 
JIM€T WHTeEPMpeTHPOBAaTb pe3yIbTaTbI IKCIICPHMCHTAJIbHbIX 
UcciIeqOBaHHH TWepeHoca HOHOB B MeCMOpaHHbIxX CHCTeMaXx, 
CCIM JaHHbIN Mpollecc MmpoTeKaeT B TasIbBaHOCTaTHYeCKOM 
pexume. Hekotopsie 3sIeKTPOKHHeETH4eCKHe MpOMecchl CBA3a- 
HbI C HOABJICHHEM PaclUIMpeHHOH OONacTH MpOcTpaHCcTBeHHOrO 
3apajja Ip CBepxnpeeIbHbIX TOKaX. OnucbIBad (POPMUpOBa- 
Hue yKa3aHHOW OOacTH, MOXKHO BbIACHHTb, KAKUM OOPa30M 
3aBHCAMe OT Hee MpOLeccbl BJIMAIOT Ha TepeCHOC HOHOB B 
TajIbBaHOCTaTHYeCKOM pexKUMe. 
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Introduction. Membrane systems form the basis for electrodialysis machines, nano- and microfluidic devices, 
which are used in water treatment, agricultural products (milk, wine, etc.) processing, performing chemical analysis, 
and in other areas of activity [1-4]. In numerous mathematical models of mass transfer processes in membrane systems 
for potentiostatic or potentiodynamic modes, the electrical mode is defined as a potential jump between two equipoten- 
tial planes parallel to the membranes. A detailed review of the papers devoted to modeling for the potentiostatic mode is 
presented in [5—7]. 

All the while, in the practice of electrodialysis, electrochemical characterization of membranes (chronopotenti- 
ometry, impedansometry, etc.), the galvanostatic mode is often used, in which a constant average current density is 
maintained at the interface. A huge amount of experimental data has been compiled about this mode. They must be in- 
terpreted [8—10]. Studies in the mathematical modeling of the galvanostatic mode are conducted in several directions. 

The first direction is the inverse problem method. As is clear from the name, this is about solving an inverse 
problem: for the specified current density at the “solution — membrane” interface, the corresponding potential jump is 
found, and then the problem for the potentiostatic mode is considered [11]. The low efficiency of this method is due to 
the fact that its implementation requires multiple solutions to the problem in potentiostatic mode for one given value of 
current density. 

The second direction is the decomposition method. In this case, the system of Nernst — Planck — Poisson equa- 
tions is replaced by a system of decomposition equations [12—16]. The assumption of a quasi-uniform charge distribu- 
tion enables to obtain a model for the galvanostatic mode in the approximation of Ohm's law [17-20]. 

The third approach can be called the direct method. In this case, the equation replacing the Poisson equation is 
derived for the current density in the desalination channel [21]. 

The galvanostatic mode can be described differently — by numerical solution to the Nernst — Planck — Poisson 
equations for an electric potential with a special boundary condition that allows the current density to be set as a param- 
eter specifying the electrical mode in the system. In [22, 23] for the one-dimensional case, the time derivative of the 
electric potential gradient was determined as an explicit function of the current density. This distinguishes the authors’ 
approaches from potentiostatic models in which the difference of potentials is set. 

This paper presents a stationary model of the ion transfer process in membrane systems for the galvanostatic 
mode. It is based on the Nernst-Planck-Poisson equation system with the boundary condition that enables to establish 
the current density as a parameter that sets the electric mode in the system. This solution is similar to [22, 23]. The dif- 
ference is that the proposed model is two-dimensional and considers the variability of the current density along the 
channel. 

Materials and Methods. The desalting channel of the electrodialysis apparatus (EDA) formed by the anion- 
exchange (AEM) and cation-exchange membranes (CEM), is taken to mean a membrane system. A binary electrolyte 
solution is pumped through it at V, average rate. 

In Fig. 1, x is a coordinate normal to the membrane surface varying from 0 (border with AEM) to h (border 


with CEM); y is a tangential coordinate to the surface of the membrane varying from 0 (channel entrance) to / (channel 
exit). 
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Exit 





0 Entrance hx 
Fig. 1. Scheme of EDA desalination channel. Concentration profiles of C,; cations (solid line) and C, anions (dot line), 


V forced flow rate are shown. 


(System of equations. Consider a two-dimensional stationary case of the system of equations describing the 
transfer of a binary electrolyte with no chemical reactions [24]: 


j S227 C.Vo—-D,VC,+CV, i=1,2, (1) 
RT 

-divj,=0, i=1,2, (2) 

E€,Ab =—F(z,C, +2,C,), (3) 

I =F(z,j, +2J,)- (4) 


Here, j , D,,z; and C, are, respectively, flux, diffusion coefficient, charge number, and molar concentration of the i- 
th ion; @ 1s electric potential; V is flow rate of the electrolyte solution; €9 is electric constant; ¢, is relative dielectric 


constant of the electrolyte solution (assumed to be constant); / is current density; F is Faraday constant; R is gas con- 
stant; 7'is absolute temperature; Ve a: T, o, C, , C, are unknown functions of x and y coordinates. 

The Nernst — Planck equations (1) describe the ion flux due to electromigration, diffusion, and convection; (2) is 
the material balance equation in the stationary case; (3) is Poisson’s equation for the electric field potential; (4) 1s cur- 
rent density in the electrolyte solution. We assume that the velocity distribution in the channel corresponds to the 
Poiseuille flow [24]: 


X X 
V=0; v, = 67,=(1-2} (5) 


To determine the unknown functions listed above, it is necessary to set boundary conditions for the system (1) 
- (3). Let us consider two electrical modes: potentiostatic, when the potential jump at the system boundaries is set con- 
stant, and galvanostatic, when the density of the current flowing through the interface is constant. 

Boundary conditions for modeling potentiostatic mode. We assume that the surfaces of ion-exchange mem- 
branes are equipotential. The system (1) - (4) includes the potential of the electric field only in the form of derivatives 


with respect to the spatial coordinate. In this case, only Ad = o(h, y)— (0, y) potential jump is significant, where AQ is 
a known function, so we set, for example: 
(0, y)=0. (6) 
Then 
Ad = const . (7) 
The conditions (6) and (7) determine the potentiostatic mode. Other boundary conditions are given below. 


At the “AEM - solution” interface (x = 0), the concentration of coions (cations) is determined from the condi- 
tion of continuity of their flow at the “membrane — solution” boundary considering the selective AEM properties [16]:. 


oF Aw _(1-5)L,0.y) 
(4 oP \o.»)- zFD,  ” 8) 


1 1 
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where 7; (i= 1, 2) are effective numbers of transfer of counterions in the membrane (CEM and AEM, respectively); 7; 
are numbers close to 1, and also for 7; = 1 ideally selective membrane, and the condition (8) turns into the condition of 


membrane impermeability for coions. 


The concentration of counterions (anions) depends on the exchange capacity of the CEM, which can be specified 
as: 


C,(0,y)=C,,, =N,Co. (9) 


Here, V, constant shows how many times this concentration differs from the concentration in the volume of the solution 
[25]. 
At the “solution - CEM” (x = h) interface for ion concentrations, the conditions similar to those at the “AOM - solution” 
(x = 0) border, are accepted: 
C.(t,y)=C,,, = N.C), (10) 
OC, F 0 1-7 )1.(A, 
. 24 20,9 (ny) =! Jk a} 
Ox RT Ox z,FD, 
At the channel entrance (vy = 0) a uniform distribution of ion concentrations is assumed: 
Ci(x0)=C,, 2=1,2, (12) 


The condition for the potential of the electric field 1s obtained from the equations (1) and (4) considering the 
absence of current flow through /, input (x, 0, ft) = 0: 





(11) 





0o(x, 0) __ 2D, OC, (x, 0) 27D: OC, (x, 0) | (13) 
Oy F(Z, D 42,D,)C, Oy Oy 
At the channel exit (vy = /), ions are freely carried by the solution flow: 
(ii, j, )= (i,-—— =D, C,Vo-D,VC,+CV)=CV,, i=1,2. (14) 


The system (14) also means that the sum of the diffusion and migration tangential components of the anions 
and cations flow is equal to 0: 





a =e" (a1) =0= 12: (15) 
Oy RT Oy 


A “soft” condition is accepted for the potential, meaning that there are no sharp changes in the potential at the 
channel exit: 


00/1) = 
(x,/) =0. (16) 


The boundary problem, which includes equations (1) - (4) and boundary conditions (6) - (16), simulates the po- 
tentiostatic mode, and its control key condition is (7). 

Boundary conditions for simulating the galvanostatic mode. Under the simulation of the galvanostatic mode, 
the condition (6) remains, and (7) should be replaced by the condition related to the specified value of i,, average cur- 


rent density at the “solution —- CEM” interface (x =h ). 
To derive such a boundary condition, we substitute the relations (1) into (4) and express the gradient of the 
electric field potential: 
RT ~ - 
Vo = —————_——_(/ + F(z D VC, +2,D,VC,)-F(z,C, +2z,C, )V). 17 
p P(eDC+2D,C)\ ( 1~1 1 2 2 5) ( 1~1 2 >) ) ( ) 
Assuming x = / in (17), we obtain the relation that connects the gradient of the electric potential with the given 
value of the current density at the boundary, that is, the boundary condition at the “solution - CEM” interface: 


[+ Fz,D, = + Fz,D, ae 
a, ©, 


& lh, y). (18) 





Ox*” F? zZDC,+2z,D,C, 


In this case, 7, current density must satisfy the condition: 


he | 
; [1..y)dy = i,,. (19) 
0 
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The mathematical model of the galvanostatic mode consists of the system of equations (1) - (4). The boundary 
conditions (18) and (19) replace the condition (7). The other boundary conditions coincide with those for the potenti- 
ostatic mode. 

The potential jump in the galvanostatic mode 1s a calculated value. 

Transformation of the boundary conditions for simulating the galvanostatic mode. The condition (19) is incon- 
venient for a numerical solution, since it contains an integral. One of this condition conversion options is given below. 

In the stationary case, the current density is a solenoidal vector. Indeed, if we multiply (2) by z,and sum up, 


then div/ = 0. Consequently, there is such 1 function that 


Map, 22. (20) 
Ox * Oy 


Using 1 function, the conditions (18) and (19) are rewritten as follows: 


: - a + Fz D, “ + Fz,D, “ 
O Vy X X 
See (fps gy ee cee a 21 
ax (1,7) F’? z, DC, + z5D,C, (sy) ae 
1 17d 1 
—[1,(h,y)dy =--[(h, y) dy = -— (nh, - 0 (h,0)) =i, (22) 
ie 1+ Oy / 
The equation (22) can be rewritten as: 
n(h, 1) — (A, 0) =-i,l . (23) 


To close the expression system, it is necessary to obtain an equation for 7 function. To this end, as in [15, 16], 


we inject a linear differential operator, which is a vortex function (rotor) in the two-dimensional case, for W arbitrary 


(W) = ( OW, } (24) 


two-dimensional vector: 





It is easy to check that: 
1) r(Vu) =0 for any u smooth function; 


2) r(uW) = (Vu, W), +ur(W) for any u smooth function and any W smooth vector. 


= O O - anes ; 
Here, (Vu,W), = Paid _“yw is skew-symmetric scalar product of Vu and W vectors, moreover, (d@,a), = 0 is for 


Oy 


any a vector. 
Applying (24) to the equation for current density (4), we obtain: 


r(f) = Far(j,)+ Far). (25) 


Using the formula of flows (1), we obtain the ratio: 
r(j,) = -<—2 DCN) —Dr(VC,)+r(CV),  i=1,2. (26) 
Hence, considering the properties of r operator: 
rj,)= -——2D (VC,,V), +(VO,V) +CrV), i=1,2. (27) 


» Ol of ; 
Considering (27) and r(/) = — — a = Ar, the equation (25) can be written as: 
X oy 


2 
An --F (a9. Sea0 22 \2_{ 2p Lazio ee), 





RT Ox * * Ox Jy éy dy } ax 
O O O O OV, OV. oe 
+F z,—b+z,—2 V -F poe V.+F (z,c,+2,c,)| —--— |]. 
Ox ce y oy Oy Ox oy 
It follows from (28) and (20) that 1 function is determined up to a constant, therefore we can assume: 
n(h,0)=0. (29) 


Then we obtain from (23) 
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n(h,/) =—i,, 1. (30) 
The conditions (29) and (30) are boundary for n function. 


Mathematical model for the galvanostatic mode in a dimensionless form. For a numerical study of boundary 
value problems, it is convenient to go to a dimensionless form. So we can simplify the equations and find out the actual 
number and set of parameters that determine the system behavior. Dimensionless variables describe a class of similar 
processes characterized by the same value of dimensionless numbers. 

Characteristic values describing the problem. When simulating mass transfer processes in the EDA desalina- 
tion chamber, a number of characteristic values are taken: 

- for spatial coordinates — fh intermembrane distance; 


- for ion concentrations — C, volume concentration of the electrolyte; 
- for speed — V, average speed of the forced flow; 
- for diffusion coefficients D =D ,D, (z,-z,)/(D, z,-D, z ,) diffusion coefficient of the electrolyte; 
- for electric potential — ~, = RT/F thermal potential; 
- for current density — i, = FDC, /h value (analogue of the limiting density of diffusion current); 
- for ion flow — j, = DC, /h diffusion flow. 
Transition formulas. We translate the equations into a dimensionless form through the following relations ( (uz) 


index denotes the dimensionless variants of quantities): 








iar, yas, =, 7 = C= i= 12, 
0 0 (31) 
go” =a. a) Se i” = 4 , j® oe ee pW Pi 
0, ly FDC, I; D 


The system of equations in a dimensionless form has the form ((u) index is omitted for simplicity): 


—_ 





j, =-z,D,C,Vo-D.VC,+PeCV, i=1,2, (32) 
—div j,=0, i=1,2, (33) 
sAo =—-(z,C, +2,C,), (34) 
An =- Ppa ae Jp. ay + 
Ox Ox } Oy Oy Oy } Ox (35) 
OV, OV, 
+Pe pag V,—Pe page V+ Pe(z,c,+2,¢,) eet eee 
Ox Of J * oy Oy Ox oy 
[=z ij, eee (36) 


The system of equations (29) - (35) contains two dimensionless numbers: the Peclet number Pe =V,h/D and 


e=6.6,RkT if (C,h°F’) . The physical meaning of ¢ parameter is that it is double square of dimensionless Debye length - 


[,: €=2(I,/h) [5]. 
Estimation of the parameter values shows that under natural conditions for electrodialysis, the Peclet number 
has 10° —10° order, ¢ number has 10°'° — 10°’ order, that is, it can be considered a small parameter. 
For computational convenience, we transform the system of equations through plugging the flux density (32) 
in the equations (33) and (36): 
div(—z,D,C, Vp—D,VC,+PeCV)=0, i=1,2, (37) 


— 


[= Sz (-z,D,C, Vo-D, VC, + PeCV). (38) 
i=] 


Thus, the system of equations contains the following unknown x, y functions: C,, C,,9, J,, Ts The fields of 
C,, C, concentrations and @ potential are determined by solution to the equations (37), (34), respectively. /,,/, current 
density components are calculated using (38). Speed distribution (5) in dimensionless form: 


V.=0, V,=6x(1-x). (39) 
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The presence of ¢ small parameter in the Poisson equation (34) means that the boundary problem is singularly 
perturbed. This significantly complicates its numerical solution, since such problems are stiff [26]. @ potential of the 
electric field and C,, C, ion concentrations change very quickly in a narrow boundary layer whose thickness is equal to 
[, Debye length [5]. 

To solve this problem, it is advisable to compact the computational grid in the boundary layer and use special 


methods for solving stiff problems [26]. 
Boundary conditions in dimensionless form. At the “AEM — solution” interface (x = 0): 

















6C ab (1-7, )1,(0,y) 
+ z,C, — |(0, y) = ———__,, 40 
( Ox a a) y) z,D, a) 
C,(0,y)=N,, (41) 
(0, v) =0, (42) 
on 
—(0,y)=0. 43 
my 65¥) (43) 
At the “solution — CEM” interface (x = 1): 
C(Ly)=N,, (44) 
ac; a (1-7, )/,.(L¥) 
+ z,C, — |(1, y) =~. 45 
( Ox 2 2) y) z,D, oo 
cules, oa D oe 
ab Gy | | ox 
—(1, y) =-—| —-——_——__—__|(1, y ), 46 
ay | y) z, DC, +25D,C, ( y) ) 
on 
—(l,y)=0. 47 
= (bY) (47) 
At the channel entrance ( y = 0): 
& (x,0) =|, 1=1,2, (48) 
ay a) D {a0 een es Jos 0), (49) 
oy Z DZD), oy oy 
n(x, 0) = 0. (50) 
At the channel exit (y = /): 
[- ag 2G. 2) s.1)=00=12 (51) 
oy oy 
(x51) = 0, (52) 
My 
1(x,/) =-i,,/, i,, =const. (53) 


After numerical calculation of the system (34), (35), (37) - (53), A@ potential jump in the desalination channel 


is determined by the formula: 
aes 
Ag=—|o(Ly)ay. (54) 
0 


The numerical solution is found by the finite element method using the Comsol Multiphysics package on an 
uneven computational grid (the density of grid elements is increased at the “solution — membrane” boundaries) [27]. 

Research Results. The calculations are performed for ¢ = 1,9-10 , Pe = 2355 , which corresponds to the fol- 
lowing values of the system parameters: 

- input concentration of the electrolyte solution of NaCl: Co = 0.1 mol/m’; 

- temperature: 7 = 298 K; 

- diffusion coefficients of cations and anions, respectively: D, = 1.33: 10° m’/s, D, = 2.05 - 10° m’/s; 

- numbers of counterions transfer in the membranes: 7, = 0.972, T, = 1; 

- ion charge numbers: z, = 1, z. =—1; 

- ratio of the counterions concentration at the boundary with the membranes to its value at the channel en- 
trance: V,=N,=1; 
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- channel width: h=10° m: 

- channel length: /= 2-107 m; 

- rate of the solution pumping: Vy) =3,8-10° m/s. 

Fig. 2 shows the fields of C; u Cy concentrations, © potential and n functions calculated at i,, = 1,57), current 
density, where i,;,. is the limiting current density determined by the Leveque formula (in dimensionless form) [28]: 


I WV.) 
Ss La me 09), (55) 





Lim 
T -t 


1 1 


Here, ¢; = 0.395 is kation transport number in the solution [9]. 
r 








Fig. 2. Concentration fields of cations C, (a) and anions C) (b), © potential (c) and 7y function (d). Model calculation for 
galvanostatic mode at 7,, =1.57,, current density. 
The forced flow acts in the channel, therefore in the areas near the membranes, the depletion of the ion concen- 
tration increases with distance from the channel entrance (along the direction tangential to the membrane surfaces). Ac- 
cordingly, the thickness of the expanded space-charge region increases along the channel (Fig. 3). 
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c) 


Fig. 3. Results of numerical calculation based on the model for galvanostatic (solid lines) and potentiostatic (dot lines) modes at 
i, = 1,5i,,, current density in cross sections y= 0,1/ (1), y=0,4/ (2), y=0,9/ (3): C, and C, concentration profiles (a); en- 
larged fragment of fig. 3, a (b); p = z,C, + z,C, space charge density (c) 

The volt-ampere characteristic (VAC, curve | in Fig. 4) is calculated on the basis of the galvanostatic model. 
In this case, the specified current density changed: i,, = 0; 0,005ijn3 ...3 1,5izm. 
40 





30 
:: 
20 
10 
0 10 20 30 AD 
Ag 
Fig. 4. VACs calculated on the models for galvanostatic (1) and potentiostatic (2) modes; ij;,, limiting current density (3) determined 
by formula (55) 


The limiting current density 1s determined by the point of intersection of the tangents to the VAC in the initial 
part and on the inclined plateau [5]. Here it coincides with i;;,, value estimated using the approximate Leveque formula 
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(55). Besides, the VAC was calculated on the basis of the potentiostatic model under a change in the potential jump Ago 
= 0; 0.4; ...; 40 (curve 2 in Fig. 4). As Fig. 4 shows, VAC 1 and 2 are congruent. Thus, there is a unique correspondence 
to each potential jump of a certain current density, and vice versa. 

Fig. 3 and 4 demonstrate a fairly good agreement of various physicochemical transport characteristics calculat- 
ed on the models for galvanostatic and potentiostatic modes. This proves the adequacy of the model of the transport in 
galvanostatic mode built by the authors. 

The computational costs of the proposed model were estimated as follows. The time spent on solving boundary 
value problems for the galvanostatic and potentiostatic modes was recorded. Herewith, the specified accuracy values of 
the calculations and the system parameters were the same, andi, =1,5i,,.. Then time indicators were compared. Thus, 


it has been found that the calculation on the galvanostatic model requires 1.6 times more time. This is due to the fact 
that: 

- galvanostatic model contains an additional equation for determining the current density distribution; 

- for the potential at x = 1boundary, the condition of the second kind (46) 1s established. 

Complex potential of the electro-membrane system for the galvanostatic mode. In [25], a generalization of the 
impedance of an electrochemical system using n function introduced under the electroneutrality conditions was pro- 


posed. A similar generalization is acceptable in this case. Note that when using electrochemical impedance, an object is 
considered only as a “black box”, and its intrinsic properties are determined indirectly. The results obtained above ena- 
ble to introduce the concept of a complex potential of the electro-membrane system: P = +i-1. Complex potential is 


a coordinate function of any point inside the object, therefore, unlike the electrochemical impedance, it allows us to 
explore the intrinsic properties of the object. 

Discussion and Conclusions. A new technique of mathematical simulation of the stationary process of mass 
transfer in the galvanostatic mode for membrane systems is described. In this case, the two-dimensional case is consid- 
ered with the use of a special boundary condition that enables to specify current density in the system. The equations for 
the current electric function are presented. The numerical solution results for potentiostatic and galvanostatic models are 
in good agreement. This shows the adequacy of the proposed transport model in the galvanostatic mode. 

The developed model makes it possible to interpret the experimental studies results of the 1on transfer in mem- 
brane systems if this process proceeds in the galvanostatic mode. Some electrokinetic processes are associated with the 
appearance of a dilated space-charge region under over-limiting currents. The description of the formation of this region 
allows us to find out how the processes dependent on it affect the ion transfer in the galvanostatic mode. 
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1,2,3 


Introduction. A special case of synthesizing the same 
electromechanical control system by the Pontryagin maximum 
principle and by the synergetic synthesis method is considered. 
The task was to solve the synthesis problem of the time 
optimal electromechanical position control system; herewith 
the travel resistance modulus linearly depended on the output 
coordinate of the system. This approach to the selection of the 
synthesis problem was because the synthesis of time optimal 
systems is one of the most widespread problems, and it is 
solved by increasing the efficiency of the existing control 
systems. 

Materials and Methods. Synthesis of the time optimal linear 
control system based on the maximum principle is a widely 
accepted problem in the modern control theory. However, the 
procedure of synergistic synthesis does not have such 
formalization. This being the case, the paper suggests an 
approach that brings together these two methods, which, in our 
opinion, will increase the efficiency of the synergistic 
synthesis method through adding some features of the 
synthesis methodology for optimal systems. 

Research Results. The paper formulates two key concepts. The 
first one is as follows: the application of the maximum 
principle for an object of the DC motor class when 
synthesizing the positioning algorithm under the conditions of 
linear loading functionally dependent on the engine rotation 
angle allows the time optimal system to be optimized. The 
second concept states that synthesis of a control system based 
on the synergistic approach enables to obtain a system close to 
optimal (quasioptimal), but after modifying the synergetic 
synthesis method itself. A hypothesis is formulated on the 
possible connection between the introduced (when 


implementing the procedure of state space extension in the 


” The research is done within the frame of the independent R&D. 


* 


ok KK v 
.PaOota BbITOJIHeHa B paMKax HHuIMaTHBHOK HUMP. 


" E-mail: i090206.lapshin@yandex.ru, tur805@mail.ru, nikaapp@rambler.ru 


JloHcKonw rocyqapcTBeHHbIM TeXHHUYeCKHM yHUuBepcuter, r. Poctos-Ha-Jlony, Poccuiickaa Denepaina 


Beedenue. PaccMoTpeH 4acTHbIit Cry4yall CHHTe3a OHO U TO 
%Ke DJIEKTPOMeXaHHyecKOH CHCTeMbI yiipaBJIeHHA MeTOOM 
MakcuMyMa IJoHTpsrMHa WM MeTOAOM cHHepreTHyeckoro 
cuHtTe3a. B kayecTBe 3a7a4un Oba onpeneneHa 3aaua 
CHHTe3a ONTHMaJIBHOU 10 ObICTPOeHCTBHIO 
9JIEKTPOMeXaHHyecKOH CHCTeMbI TMO3HIMOHHpOBaHHA, IIpu 
9TOM MOMECHT COMPOTHBJICHHA JIBHXKCHHIO JIMHeCMHO 3aBHCell 
OT BbIXOHOM KOOPAMHAThI CHCTeMBI. DTOT MOAXOA K BbIOOpy 
3aqa4uH CHHTe3a Obit OOyCcOBIeH TeM, 4TO CHHTEe3 
ONTHMAJIbHBIX 10 ObICTPOeCHCTBUIO CHCTeM ABJIAeCTCA OHOU 
M3 CaMBbIX INMpOKO paciipocTpaHeHHbIx 3aa4, KOTOpaA 
pelaeTca IIpw MOBbIMeHHH 9d@eKTHBHOCTH JeMCTBYFOIINX 
CHCTeM ylipaBJieHHa. 
Mamepuanot u  memodvi. CHHTe3 ONTHMabHOK m0 
ObICTPOJeCHCTBHIO JIMHEHHOU CHCTeMbI yiIpaBsIeHHA Ha OCHOBe 
IpwHuwuna MaKCHUMyMa — IIMpOKO paciipocTpaHeHHas 3ayaua 
B COBPeMeHHOM TeopHu ypaBlIeHus. OHako mpollesyypa 
CMHepreTHyuecKoro CHHTe3a TakOM PopMasM3allun He UMeECT. 
Ucxoja 43 STOFO, B CTaTbe TIpeIOHKeH NOAXON, cOmMKaIOMIM 
9TH TBA MeTOa, KOTOPBIM, 10 MHeHHIO aBTOPOB, MO3BOJINT 
NOBbICHTb 3(@eCKTHBHOCTbh MeTOa cCHHepreTHyuecKoro 
CHHTe3a, 2O0aBHB B Hero HeKOTOpbIe OCOOeHHOCTU 
MeTOJOJIO‘HN CHHTe3A OIITHMAJIbHBIX CHCTeM. 

Pe3zyiobmamoi ucciedoeanua. B padote cPbopMyIMpOBaHBbI Ba 
OCHOBHbIX Hay4HbIX WoOKeHHA. Ileppoe - mpumMeHeHue 
IpvHuuna MaKCMMyMa Jif OObBeKTa kuIacca JIBHTaTesib 
MOCTOAHHOTO TOKa lIpu CHHTe3e ayIrOpuTMa 
NO3SHIMOHMpOBaHHaA B YCJIOBHAX JIMHeEMHOM Harpy3Ku, 
(YHKIMOHAIbHO 3aBNCAeH OT yrula MOBOpOTa ABuraTeJ4A, 
MO3BOJIACT ONTHMH3UpOBaTb CHCTeMy 110 ObICTpOAeHCTBUIO. 
Bropoe - CHHTe3 CHCTeMbI ylipaBsIeHHa Ha OCHOBe 
CHHepreTH4uecKOrO TOAXONa MO3BOJIAeCT MOUYAHTb CHCTeMy, 
ONM3KYIO K ONTHMAaIbHOM (KBa3HONTHMasIbHYy!O), HO yKe 
mocue MOWM@uKaNMH CcaMOro MeTOJJa cCHHepreTHyueckoro 
cuHTe3a. CiopmMysMpoBaHa rumoTe3a O BO3MO2%KHOM CBA3H 
IIpoueypbl 
COCTOAHHUA «60S sBB-—SC Mee 


MOXY 
paciumpeHnna 


BBOJUMBIMH, pu  peasM3aluu 
IIpOCTpaHcTBa 


CHHepremTuvyecKOrO CHHTC3a, WOCTOAHHbIMU BPCMeCHH C 
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synergetic synthesis method) time constants with the optimal 
switching time of control defined in the maximum method. 


Discussion and Conclusions. The synthesis through the 
maximum control technique and the ADAR method is 
performed. In virtue of the comparison of efficiency of these 
methods, a hypothesis is put forward on the possible 
compatibility of the studied methods. 


Keywords: maximum principle, optimal control, operation 
speed, control algorithm, synergetic synthesis, ADAR method. 
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Introduction. In the 20th century, the engineering requirements, in particular, in space engineering, put 
forward a range of problems for which a new theory was developed — the theory of optimal control [1]. One of the main 
techniques for the synthesis of optimal control systems is the maximum principle developed by Soviet mathematician L. 
S. Pontryagin and his disciples in the fifties-sixties of the 20th century [2]. The application of this principle is based on 
the formalization of the synthesis problem with the transition to the form of the Mayer problem and the subsequent 
solution to systems of linear or, 1n some special cases, nonlinear differential equations [1]. It is worth noting that the 
need to solve systems of differential equations, and in the nonlinear case this is not always possible, is in many ways a 
limitation of both the maximum principle itself and the whole concept of synthesis of optimal control systems. 

In recent decades, a new approach to the synthesis of control systems has been widely adopted. It relies on the 
synergetic concept of the analysis and synthesis of systems. The technique used, the author of which is A. A. 
Kolesnikov, is called the method of analytical design of aggregated regulators (ADAR) [4-5]. The proposed approach is 
based on the concept of synthesis of nonlinear feedbacks. They provide the asymptotic stability of the control system 
with respect to the required motion of the attractor in the state space [4]. This method differs from the methods of 
synthesis of optimal control systems in the absence of both the optimization criterion of the control system and the 
statement of the synthesized control optimality. From the point of view of the implementation of the synthesis 
procedure, the ADAR method has an undeniable advantage over the synthesis methods of optimal systems, which is 
expressed in the absence of restrictions on the nonlinearity of the system of differential equations [4-5]. 

From a practical standpoint, synthesis of the process or object management should be able to answer the 
question if there is another control that has the property of superiority over all others. Thus, the modern mass production 
constantly requires efficiency improvement to ensure market competition. One of the most popular ways to increase this 
efficiency is to optimize management processes. This approach is applicable both to the systems of automated assembly 
of equipment [6—7] and to the metal-cutting systems, in particular, to drilling control systems in metal-cutting machines 
[8—15]. The idea of combining these approaches to synthesis of the control systems has, from the authors’ the point of 
view, an undoubted practical effect. In science terms, it 1s important to combine the ADAR method advantages, which 
are expressed in the possibility of considering synthesis of the complex nonlinear process dynamics in mechanical 
engineering [16—19], with a neat and definite formalization of the synthesis problem formulation and assessment of its 
achievement in the maximum principle [1—2] . 


1. Synthesis of basic mathematical model and formulation of research problem 

In the modern economy, the direction of the time optimal system synthesis, which enables to obtain a 
significant increase in the global system efficiency [7—8], has become a frequent practice. Thus, under solving the 
problems of automating the assembly processes of various equipment, the task of attaching different types of parts to 
each other [6—7] often arises. Here, the economic efficiency of the entire production process depends on the speed of 
this operation. The same situation is observed in metal cutting systems on metal-cutting machines [9-19], in which the 
faster the machining process is, the lower the costs of the entire production process. Based on the reasoning, under the 
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assumption that it is necessary to synthesize a time optimal control system, which, considering the similarity of the 
processes of automated fastening of parts and drilling deep holes, can be illustrated by the following diagram (Fig. 1). 


C 


DC motor 





Fig. 1. Illustration diagram 


Fig. 1 shows a DC motor performing either an equipment fastening operation or a feed in a drilling unit. In the 
first case, C is the combined load-deflection characteristic of the material that 1s compressed under twisting. In the 
second case, C characterizes a linear increase in the feed resistance under the accumulation of chips in the chip removal 
channels of the drill. 

The electric drive, which ensures the conversion of electrical energy to mechanical energy of motion of the 
system actuators, is the basis for the support subsystem of the control system in both cases considered. Let us assume 
that in the present case, such a transformation is carried out by a DC motor with collector control, which is described 
through the following system of equations [20]: 

U-—c,0=L y + Ri 

dt 

clad 2 +M. 
dt 

where U is voltage applied to the engine manifold; i is current consumed by the motor; R, L are parameters of the 
electrical part of the engine; / is parameter characterizing the inertial properties of the engine rotor, the reduced 
inertial moment of all rotating masses; w is engine rotor speed; M,is external moment of resistance; Cy, C. are 
mechanical and electric engine constants. In this case, the moment is presented as a function of the angle of rotation of 


the engine rotor, that is M@, = Ca. With this in mind (1), we rewrite as: 


et ae 
dt 








d 
BJ aCe. (2) 
dt 
do 
—_ = 
dt 
We express the current value in the second equation, apply the obtained value in the first equation, and get: 
l LJ d°o RJ do LC RC 
U—-a=——; + —+——wo+——a 
C, Cea €c. at Ce. Ce, . (3) 
do 
—}$ =f) 
dt 
We solve the first equation with respect to the highest derivative and rewrite it with the second one, and then 
we receive: 
do 
ds 
d°o Rdwm cc, C RC ae (4) 
eS 0 0 
dt Os LI J LJ LJ 
do ; R cc, C RC C. 
Let us puta =x,,@=x,,—=x,; as constants, we introduce — = a,,,(“*++—) = 4a,,,— =4a,,,— =D. 
dt L LI J LJ LJ 


Then the system (4) takes the form: 
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dx, 
— =X, 
Va 
X5 
di 3 (5) 
X3 
— = —,,X, —Ay,X, —Ay.X, +bU 
dt 31°"1 32°"2 33°"3 
or in matrix-vector form: 
x= Ax+Bu, (6) 
0 1 0 (0) Da (0) 
where A=| 0 0 L ig BH Ole HH Xo 1S 0 
—a —a —a b x U 


31 32 33 
We take the DC motor whose parameters provide the following constant values: a3,;=4.65, a32=4.6, a33=2, as a 


drive ensuring the system motion. 


0 1 0 
A matrix, with this in mind, will take the form: 4=| 0 0 1 |. From knowing A matrix, we find A 
4,65 -4,6 -—2 


eigenvalues of A matrix: A, = —1.2666 + 0.00001, A, =—0.3667 + 1.88061, A; = —0.3667 — 1.88061. 

As is clear from the obtained eigenvalues of the matrix, the control system is asymptotically stable according 
to Lyapunov [21]. Moreover, we can argue about the oscillatory nature of the processes proceeding in the system, since 
the eigenvalues contain not only negative real parts, but nonzero imaginary parts. 

2. Synthesis of control by Pontryagin’s maximum principle 
First, we formulate the problem of the optimal control synthesis in the following form: 

- using the maximum principle for the object described by the System (5), determine the optimal equation 
algorithm that ensures the transfer of the object from the initial state x,(0)=50, x.(0)=0, x3(0)=0 to the final state 
x}(T)=0, x2(T)=0, x3(T)=0 for T minimum time. Herewith, it is necessary to determine transition count, switching 
torque, and to construct curves of u(t) control and x;(t), x(t), x3(t) coordinates. 

Jul\<Umax Yestriction is imposed on the control action. The system parameters are as follows: 
a,, =1,a,, =12,a,, =20,b=1, Unax=440 V. 


Solution: We write the Hamiltonian: 
2 
i= wl ? (7) 
i=l 


where w, and f; are coordinates of w and 7 vectors. In addition to the system (5), we construct a system of equations 


for wy, auxiliary variables, where i = 1,2 , using the following relation: 


dw, 2. Of 
Vi =-)'y, Ji 
jal O% 








r= 17. 
dt O 


i 


Or open : 


dW, of OP Ss 
st — |] tt Ae fe Se 
dt 2 ™ ax, me ax, - 





dW, oh oS os 
= =-| —y, +——y,+— 8 
dt E 7 ax, a ax, ” ®) 
d 0 0 0 
Ue i oy iy, 
dt Ox, Ox, Ox; 
Considering (5), the equation system (8) takes the form: 
aaa =4,65¥, 
dt 
Y 
— =-Y +4,6¥, (9) 
i 
cco —P, +2, 


dt 
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Y 
cas =4,65¥, 
dt 
d'¥ 
—s =P +4,6¥, (9) 
Y 
— -=-P, +2, 
t 
Eigenvalue matrix for (9) case: 
1.2666 0 0 
Ar? =| 0 0.3667 + 1.88067 0 
0 0 0.3667 + 1.8806i 
The solution to the system (9) for the diagonalized case of A¥~? matrix takes the form: 
WY = Ce! 2666 
1 “1 
w _ C ge ON ree) (10) 
| es: : 


We = Ce 03667-1.88061)1 
3 3 


Considering V” matrix, the solution in the initial basis will be: 
Y= —0.9474C,e°" + (0.5812 —- 0.334 71)C, ero E80 + (0.5812 + 0.33471) Cee 
Se | 1 893C,e'°°"" ie 0.6883C, e007 S800! = 0.68830, e077 18800! 
WY, =—0.2581C,e°°" + (0.1812 + 0.20861)C, e070 "" + (0.1812 —0.20867)C,e0 


where C,, Cz, C3 are integration constants. 
The general expression describing the Hamilton function: 
H=P f+ Poi +f (11) 
In the expression (11), an important — from the point of view of the method of synthesis — role is played by the 
member, which includes the control: 
H’ =,U =(-C,0.258 le! + C, (0.1812 + 0.20861) + C, (0.1812 —0.20861)E rrr "BU" 
In order for H Hamiltonian calculated by the formula (11) to take the maximum positive value, H term must 
be always positive and the greatest. For this, the optimal control algorithm should be u(t)=o0U,,,,, where: 
o = sign(—0.2581C,e'°" + (0.1812 + 0.20861) E07" +. (0.1812 — 0.20861) Err TO) 
u(t) optimal control is a piecewise constant function, taking +U,,,, values, and it has no more than two 
intervals of constancy, since the nonlinear function 
—0.2581C,e'°°" + (0.1812 + 0.20861) E07" +. (0.1812 — 0.20861) E77 18°" 
changes the sign no more than once. In this case, the possible sign change occurs from plus to minus, that is, to fulfill 
the maximum principle, it is required to first apply U= +U,,,,, to the engine, and then U= —U,,,,. Let us verify these 
arguments through constructing ‘¥, obtained functional dependence in the Matlab package (Fig. 2). 





beer O.1 O.2 O.5 O.4 0.5 0.6 OF 0.0 O.9 1 


Fig. 2. Graph of ‘Y, variance 
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As can be seen from Fig. 2, ‘¥, does change the sign from (+) to (—) only once. With this in mind, we define an 


optimal equation algorithm that ensures the transfer of an object from x,(0)=0, x,(0)=0, x3(0)=0 initial state to 
X,(T)=50, x2(T)=0, x3(T)=0 final state for 7 minimum time. Here we note that under modern conditions, there is no 
need to obtain an analytical solution to the original system of equations. Using the available capabilities of modern 
software packages like Matlab, we straightforward and clearly can obtain a numerical solution to the case under 
consideration. 

The numerical simulation results for the calculated time optimal control algorithm are presented in Fig. 3. The 
model parameters are selected in such a way that it fulfills the boundary condition required by the right-hand border. 

Xi, rad 


ol) 
40 





0 0.5 1 15 & See 





a5 1 15 
b) 


Fig. 3. Graphs of output coordinate (a) and control switching 
at t, = 0.5 (b) 
Fig. 3 shows that the control switching time from (+U max) to (-Umax) is selected as ¢, = 0,5 seconds, and the 


total control time is T=1.3 seconds. In this case, the control system comes to the required output level, that is, the 
right-hand boundary of the boundary conditions is reached. 
3. Synthesis of control system through ADAR method 
To synthesize a control system using the ADAR method, it is possible not to make a transition to the abstract 
case of the state space (see (5)), but it is easy to use the original system (2). However, it needs to be translated to the 
following form: 


dx, 
i 
wy) 
— =—-@,,X, + a,,x : (12) 
gt 21°71 207.8 
X. 
— = —d,,.X, — dy,xX, +bU 
dt 
, C, C c R 1 
where Q=x,,@=xX,,i=x,are variables; ra. Ay; - ie = Ay ae Ay, ae b are introduced as _ constants. 


Considering the previously defined values, we receive: a,, =1.515,a,, =2.325,a,. =1.5,a,, =2,b=1 and the system 


(12) takes the form: 


dx _ 

a 

M2 = -2,325x, +1.515x, , (13) 
t 

d 

ae =—154, 2 255 EU 

dt 


With this in mind, 4A matrix takes the form: 
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0 | 0 
Az=|-2.325 0 - 1.515]. 
0 -15 -2 
From knowing A matrix, we find i eigenvalues of A matrix: 11 = —1.2674 + 0.00001; A2 = —0.3663 + 1.88011; 
X43 = —0.3663 - 1.88011. As is clear from these results, it is referred to the control system case considered in the 


synthesis through the maximum principle. 
The coordinate characterizing the angle of rotor rotation (x) is the output coordinate of the system. Therefore, 
in order to form requirements for the desired system behavior in the state space, we introduce a macrovariable of 


Y =x, —X, = Oorder, where x,, is the specified target value of x, coordinate. The velocity value of the angle change 
must ensure that the angle of rotor rotation tends to the value specified above. To this end, we introduce another 
macrovariable — ‘VY =x,—0,(x,) => 0, where o,(x,) is a certain function describing x, > ,(x,) tendency in the 
stationary state. Further, we introduce a macrovariable of order in VY =x,—6,(%,,x,)=>0 coordinate, where 
;(x,,x,) 1S a certain function that describes the relationship between the coordinates in the stationary state of the 
system and, therefore x, > 0,(x,,X,). 

For newly introduced 'Y,, ‘¥,, ‘¥, macrovariables, we require the asymptotically stable law of change, that is, 


execution of the following system of differential equations: 











d'¥ 
T —+¥,=0 
t 
d'¥ 
T, —2+¥, =0, (14) 
t 
d'¥ 
T, —++¥,=0 
t 


where 7, T7,, T, are time constants that determine the vanishing rate of the introduced macrovariables, or, in other 


words, the shrinking rate of an arbitrary path of the system in the state space to the desired attractor. 

The system of equations introduced in the expression (14) has expanded the system state space from the 3rd to the 
6th, which is one of the main provisions of the ADAR method. Further synthesis of the control system is reduced to a 
stage-by-stage process of decomposition (compression) of this space to the initial level [1—2]. We start the 
decomposition process with the substitution of the system (14) to the first equation instead of ‘, macrovariable of its 
value in the coordinates of the controlled process. As a result, we obtain: 


dx 
Re tix) = 9, (15) 
Considering (13), the equation (15) is written as: 
Tx, +(X, -Xq,) = 0. (16) 


Considering the next step of the synthesis algorithm, because of which it was assumed thatx => 0,(x,), we find 


, (x, ) value as: 


a (x, = X1) 
(x)= r (17) 


1 
With this context, ‘¥, macrovariable is as follows: 


Veer (18) 


1 
Validity of the expression (17), in terms of the control objectives, is confirmed by the fact that the value in the steady 


state is as follows: x, > 0,(x,) => 0. Then, considering (18), the second equation of the system (14) will take the 


following form: 
dx do, (x, ) 
a ee (19) 
Or, considering 0, (x, ): 
dx, te fT, dx, a ts (=o) —(Q. 


20 
aT i ar (20) 
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At the next step of the synthesis algorithm, using the same reasoning as before, and the fact that we 
determine x, > 0,(x,,x,) from (20), considering (13), ,(x,,x,) value in coordinates of the controlled process: 


lL, — 
tte oes 066 2 066 2 , (21) 


1 Ea 
Having obtained ,(x,,x,) value in the coordinates of the controlled process, we can determine the value of 


’; macrovariable in the coordinates of the system state and solve the third equation of the system (14). 


~ Eh x, dm) =0. (22) 


After substituting the previously obtained 9,(x,,x,) values from (21) into (22) and using the system (13) at 


dx 
5a aes 
al? 


this step of decomposing the state space of the control system, we define the control value in the coordinates of the 
controlled process: 


GQht+hG +E) GBt+hB+hh) 
(ee © CLL, 
(T,+T,+T,)_ 


Ue 2 og 
Hee ee 


The expression (23) determines the asymptotically stable control action on the control system described by 
the equations (13). 


1,52, 42x, 1, 55994, 1552, 
(23) 


The simulation results of the system of differential equations (13) with the required value of the engine rotor 
angle: x;9=50 rad, and the values of 7, =0.1,7, =0.2, T, =0.3 time constants entered are shown in Fig. 4. 
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Fig. 4. Simulation results with control at 7, = 0.1, 7, =0.2 , 7, = 0.3 : transition process on output coordinate (a); 


control change schedule (b) 


Fig. 4 shows that, in comparison with Fig. 3, the control structure is the same in both cases, that is, the point 
is that at the beginning, a positive control is applied, and then its sign changes to the opposite one. The settling time 
increases dramatically to the value close to 2.5 seconds, but with this, the maximum control reaches values greater 
than 5,000 V. This is unacceptable according to the maximum method. To limit the maximum control value, we 
introduce a nonlinear link of the “saturation” type with the same threshold as in the case of the maximum into the 
control loop. The results of the system simulation, with this threshold and with selected values of 7, T7,, 7, 


constants, are presented in Fig. 5. 
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Fig. 5. Simulation results with control at 7, = 0.15, 7, =0.16,7, =0.17: 


output coordinate transition (a); control change schedule (b) 


Fig. 5 shows that the introduction of a nonlinear constraint on the control of the “saturation” type enables to obtain the 
desired result from the point of view of the control admissibility. In this case, the control structure determined by the 
maximum method is really observed, but the control is not time optimal. It should be noted here that such a 
modification of the ADAR method could also lead to the loss in system stability (Fig. 6). 


Xi, rad 


L sec 


t, sec 





Fig. 6. Simulation results with control at 7, = 0.015 , 7, =0.016,7, =0.017: 


output coordinate transition (a); control change schedule (b) 


However, despite the result presented in Fig. 6, the addition of amplitude constraint of the control signals to the 
control system obtained by the ADAR method essentially brings it closer to the solutions obtained through 
implementing the synthesis procedure by the maximum method. 

Discussion and Conclusions. Thus, for the basic model of the control system (2), the time optimal control was 
synthesized by the maximum method. For the same case, in the third part of the paper, with the same system parameters 
and boundary conditions, the control of the ADAR method was synthesized. Both of these controls look different. 
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In case of the maximum method, this is a software piecewise-constant control where the optimal time instant of 
the sign change of the control signal 1s important. In ADAR case, it is a continuous smooth functional control 
dependence on the coordinates of the system state. However, their implementation shows the structural similarity 
expressed in changing the control sign. After the ADAR method modification, the structural similarity has become even 
more visible. Note that the parametric similarity of the implemented controls can be achieved through selecting the time 
constants introduced by the ADAR method during the synthesis. This allows the authors to formulate the following 
hypothesis. 

- Selection of the values of 7,,7,,7, time constants introduced by the ADAR method in the synthesis 


procedure, in case of modification of the obtained control by the method proposed in the article, makes it possible to get 
time optimality of the synergetic control. The authors failed to obtain a general proof from this intuitive hypothesis. 
However, this task was not set within the framework of this paper. The numerical experiments with 7,,7,,7, variation 


helped us to achieve quasioptimality of the control obtained by the ADAR method. 
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HapaseubHOe NocTpocnHe TBOMUHOrO epeBa Ha OCHOBE COPTHPOBKH 
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Pocculickaa Deeparna 


Introduction. Algorithms for the parallel binary tree construc- 
tion are developed. The algorithms are based on sorting and 


N 


described in a constructive form. For the element set, the 


T(R)=O0) ang TR) = Oleg N) 


time complexity has esti- 


2 
mates, where R=(N-N)if2 


is the number of processors. 
The tree is built with the uniqueness property. The algorithms 
are invariant with respect to the input sequence type. The work 
objective is to develop and study ways of accelerating the 
process of organizing and transforming the tree-like data struc- 
tures on the basis of the stable maximum parallel sorting algo- 
rithms for their application to the basic operations of infor- 
mation retrieval on databases. 

Materials and Methods. A one-to-one relation between the 
input element set and the binary tree built for it is established 
using a stable address sorting. The sorting provides maximum 
concurrency, and, in an operator form, establishes a one-to-one 
mapping of input and output indices. On this basis, methods 
for the mutual transformation of the binary data structures are 
being developed. 

Research Results. An efficient parallel algorithm for construct- 


ing a binary tree based on the address sorting with time com- 


plexity of EO) is obtained. From the well-known 
analogues, the algorithm differs in structure and logarithmic 


estimation of time complexity, which makes it possible to 


achieve the acceleration of ON ye a21 


order analogues. As 
an advanced version, an algorithm modification, which pro- 
vides the maximum parallel construction of the binary tree 
based on a stable address sorting and a priori calculation of the 


stored subtree root indices is suggested. The algorithm differs 


in structure and estimation of 7“) =? time complexity. A 


similar estimate is achieved in a sequential version of the mod- 


* PaOota BBINOHeHa B paMKax HHUWMaTuBHO HUMP. 
™ E-mail: romm@list.ru, denchabanyuk@gmail.com 
“™ The research is done within the frame of independent R&D. 


Beedenue. Pa3spadboTaHbi asIrOpHTMbI MapasWiesIbHOrO MOCTpo- 
CHUA JBOMYHOrO JepeBa. AJITOPHTMBI BbIMIOJIHCHbI Ha OCHOBE 
COPTHPOBKH UH OMMCaHbI B KOHCTPyKTHBHOHM Popme. Ja MHO- 


N 


XKCCTBAa U3 IICMCHTOB BPCMCHHadA CJIOKHOCTbh UMCCT OLCH- 


T(R) =O), T(R) = O(log, N) 


KH , THe YACIO IWpowecco- 


2 
pos “ Se JlepeBo CTpoutca CO CBOMCTBOM eJMH- 
CTBeHHOCTH. ANTOPHTMbI HHBapMaHTHbI OTHOCHTeJIbHO BUA 
BXOJHOM mMoceqOBaTenbHOcTH. Llembro padoTbI sABIIATIACb 
pa3paO0oTka u UcclIeqOBaHHve CMOCOOOB ycKopeHusA mpolecca 
OpraHn3alun UW WpeoOpa30BaHui JIPeBOBUAHbIX CTPyKTyp 
JJAHHBIX Ha OCHOBEe aJITOPHTMOB YCTOMYHMBOM MaKCHMaJIbHO 
TlapaswIeIbHOH COPTHpOBKH JIt UX TpWYMeHeHuA K Oa30BbIM 
ollepalluamM HHDOpMalMOHHOLO MONCKa B 0a3ax JIAHHBIX. 
Mamepuaibi u memoodvi. B3au4MHO OJJHO3HAYHOe COOTBET- 
CTBH€ MHO2KeCTBAa BXOJIHBIX IJICEMCHTOB HW MOCTPOeCHHOTO JIA 
Hero J{BOMYHOrO epeBa ycTaHaBIMBaeTcA pH MOMOMM 
YCTOHYMBOH alpecHow copTupoBKH. CopTupoBka oOsIazaeT 
MaKCHMA@JIbHbIM [apasIesIM3MOM, B OlepaTOpHOH dopme 
yCTaHaBJIMBaeT B3AaMMHO OJHO3HAYHOe COOTBeETCTBHe BXOJI- 
HbIX WM BbIXOJHbIX HHAeKcoB. Ha 3sTOM OCHOBe pa3pa0aTHIBAa- 
FOTCA MCTObI B3AMMHOTO TipeoOpa30BaHHA JIROMYHbIX CTPyK- 
TYp JIAHHBIX. 

Pe3yibmamol ucciedoeanus. Wonysen 9d@eKTHBHBIN Mapas- 
JICJIBHBIM aITOPUTM MOCTpOeHHA JBOMYHOrO JepeBa Ha OCHOBe 
ayIpeCcHOou CJIO2KHOCTBEO 


COPTHpOBKH C  BpeMeHHOW 


T(N )= O(log, N). OT W3BeCTHBIX aHaJIOrOB aJITOpHTM OTJIN- 
yaeTCA CIPyKTypou HW WOrapudMuyeCKON OLLEHKOM BPeEMeCHHOU 
CJIOHKHOCTH, MO3BOIAFOIeM MOCTHTaTb yCKOpeHHA aHasIOroB 
nopsyka O(N’), a2>1. B KauyecTBe ycoBepilieHCTBOBaHHOTO 
BapHaHTa IipesoxeHa MOAMMUKalluA alropuTMa, oOecreun- 
Barlolllad MAaKCMMAJIBHO MapasiesIbHOe MOCTpoeHve JBOMYHOTO 
yjepeBa Ha OCHOBeE YyCTOMYMBOHM afpecHow COpTHpoBKH U 
alIpHOpHOrO BEIYMCJICHHA XpaHUMBIX MHACKCOB KOpHeH MOI- 
jjepeBbes. AsIrOpuTM OTIMYaeTCA CTPyKTypoH UU OleHKON 
BpeMeHHOHM cioxKHOCTH 7(1)=O(1). AHanorM4Had oOleHKa 


JJOCTMTaeTcA B MOCNeOBaTesIbHOM BapHaHTe MOMUIMpO- 
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ified algorithm, which allows obtaining the acceleration of 


O(N’) a>] 


known analogs of order. 


Discussion and Conclusions. The results obtained are focused 
on the creation of effective methods for the dynamic database 
processing. The proposed methods and algorithms can form an 
algorithmic basis for an advanced deterministic search on the 
relational databases and information systems. 


Keywords: data structures, data processing algorithms, binary 
tree, algorithms for parallel sorting. 


For citation: Ya.E.Romm, D.A. Chabanyuk. Parallel 
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BaHHOrO asIrOpHTMa, YTO MO3BOIAeT JOCTHTaTbh YCKOpeHuA 
W3BeCTHbIX aHaIOrOB HopsyKa O(N’), a>1. 

O6cyocoeHua u 3akiioyeHue. Ilomny4eHHbie pe3yiIbTaTbl 
HanpaBJICHbl Ha OpraHu3al{u1o IPeKTUBHBIX CHOCOOOB JMHa- 
MHYeCKON OOpadoTKH 6a3 WaHHbIx. IIpensox*xKeHHbIe CMOCOOBI 
WM asITOPHTMBI MOFyT COCTABHTb asIrOPHTMMYECKYI0 OCHOBY JIA 
YCKOpeHHOrO JJeTePMHHUPOBAaHHOrTO TOUCKa B PeJIAIMOHHbIX 


0a3ax JJaHHBIX 1 MH(OPMAallMOHHBIX CHCTeMaXx. 


Ku1104ueBble CJIOBa: CTPyKTYpbl aHHbIX, asIrOpHTMbI OOpa- 
OOTKH JaHHBIX, JBOWMYHOe epeBo, asIrOpuTMBI MapaliiesIbHOU 
COPTHpOBKH. 


Oopazeu Osa YUUumupoeanua: Pom, A. E. IapasnembHoe 
TMOcTpoeHve JBOMYHOrO jepeBa Ha OCHOBe COpTHpoBKH / 
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Introduction. There is a tendency to the convergence of parallel information processing technologies and various 
processor architectures in the field of modern high-performance computing. Despite the variety of processor architectures and 
ways of presenting information, the idea of parallel processing is one of the most important tasks of computer science to 1n- 
crease the data-rate. To accelerate processing speed, the authors propose to use a stable address sorting algorithm with maxi- 
mum parallelism. 


Method of parallel construction of a binary tree. For 4 = ( D, , Disses a,.) array, the comparison matrix is developed 


according to [1,2]. a,, element of this matrix is defined as 


ty. 2 Sd, 

a,, =sign (a,—a,)= 0, a, =a, 

-, a, <a. 

J i 

2 
where 7, j=1, 2, ..., 7. 
a. element in C=(c so. ea) sorted array gets the number k=) \a,,, where a,,>0 ati<j,a,,>0 at 
J 0 1 n-1 = ij ij ij 


i> 7. All comparisons are mutually independent; the sorting is stable and as parallel as possible with the estimate of 


Zz 


Suppose we are given a set of N elements, all elements of which are represented as a single-dimension array. 


N’-N ; , , 
T (==) = O(1) time complexity. On this basis, you can perform a parallel construction of a binary tree [3, 4]. 


On the set, < ordering relation is assumed. It is required to convert the array into a binary tree. For this, the described 


array sorting is performed. C’ medial array cell has 7, = i] index and is taken as the root of the tree [3]. All C array 


components to the left of C’, form a left subtree (left subarray). The components to the right of C, form a right subtree 


(right subarray), The left subarray is_ interpreted as a new = array. It similarly contains 


1(| N mee 
Feet (4 1) = a root index. Here, C, ____ is the left-nearest descendant of the root of C, tree. All 


components of the subarray to the leftof C, | donotexceedC, __ ; all components of the subarray on the right are 
N|] /1([N ee ee oe 
not less than C, . Simultaneously, the root index of 7, as.v2 =| |+} =] | = |-!]|=/7, +| —— | tight subarray 
ep. nes. 1/2 cp. pas. 9) a 2 cp f) 
is determined. At this, C, _ is the nearest descendant of the root of C, tree. The process recursively resumes in 


each pair of the adjacent subarrays: 


Romm Ya. E. and the others. Parallel construction of binary tree based on sorting 


Pomn A. E. u op. Tapasmiensnoe nocmpoenue 06ouunozo Oepeéa Ha OCHOGe COPTHPOBKU 
eee ——————————————————————— SSS SS SSS SSS aaa! 


- _ Divas ye > | 
[en v2i,1 2 ? 


. os Ie ae i 7 I 
J c 22° J V2 = ? 
p. eB. , cp. Jes. p) 


7 os: Mas eae yor _ i" yor? Zz | 
Tce 1/2',1 i io 2 : 


J. i-l = 4 i-2 —] 
Deum an ne yo a | 7 ’ i=] ’ 2 ‘wits log, N. 


As a result, all components of the lower level of the binary tree are formed in O(1) time. The process can con- 


tinue until log, NV exhaustion of the levels of the binary tree. 

The number of algorithm steps for constructing the binary tree in a parallel form is the sum of the sorting step 
and the step sequence when calculating the indices of the roots of subtrees. From here, 
T(R)=log, Nt+ t= O(log, N), where R is the number of processor elements, tis the time of binary comparison, 
and Tis the time for calculating one root index. R number of processors is determined by the maximum JN parallel 
sorting of input elements, and then — through the calculation of indices with doubling by the number of tree levels. 


When calculating the indices, this number will not exceed 2°’”" = N/2, therefore the number of processors involved 


2 





in sorting is sufficient. As a result, R will be less than [3]. Finally, the time complexity of the parallel algo- 


rithm for constructing the binary tree will be 





[= = O(log, N). 


Example [3]. The binary tree for an array of 15 elements 
xX = (14, 9, 24, 7, 11, 20, 28, 3, 8, 10, 13, 17, 21, 25, 30) is constructed as follows. 
The result of the sort is the array 


C =(3, 7, 8, 9, 10, 11, 13, 14, 17, 20, 21, 24, 25, 28, 30) 


, . , 1 
The root of the binary tree is the medial element of C array: j,, -|5|-8. C,=14 The left subarray has 


—| as 
i nee. 1/2 = Stl root, C, =9 element is the root of the left subtree, which is the left-nearest descendant of Ci. 


=] ; 
medial component. The right subarray has Ie page a+| =12 root, C, =24 element is the root of the right 


; 4-] , 
subtree and the right-nearest descendant of C, root. Further, 7 = is ] = 2, C, =7 element is the root of 
ep cp. mes. 1/4,1 5 oy 


the subtree on the left and the left-nearest descendant of the root of C, ___ subtree. In the right subtree, the root has 


Aad . 
J. V4.2 =4 | = 6 number, C, =11 element is the root of the right subtree and right-nearest child of C, 
P. acc : i) Joep. nes. 1/2 
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pan.1/2 Top. mpas. 1/4,1 


12-8-1 , 
Similarly, to the left of C. = 12-28-10 root is determined, C,, = 20 element is the root of 45] 
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the left subtree from it and the left-nearest descendant of C, _ subtree root. For the subarray, adjacent to the right one 


12-8-1 
discussed above, the root has : ae ia 124] So a number, C,,=28 element is C, __ right-nearest 


descendant and the right subtree root. The lower level of the tree will be formed by the descendants remaining to the left 


and to the right of each of four identified roots (Fig. 1): 


(0 level 






Joep. nes. 1/2 — 4 Joep. mpas. 1/2 =12 


1 level 





=14 





Jon. mpas. 1/4,2 


=a 





Jon. mes. 1/4,1 J op. mes. 1/4,2 = J op, mpas. 1/4,1 — 












Joep. mes. 1/8,1 =i J ep. nes. 1/8,2 3 Jp. mes. 1/8,3 =5 J op. mes. 1/8, 4 = a pas. 1/8,1 = Tes mpas. 1/8, 2 =11 Sx mpas. 1/8, 35 Fep. mpas. 1/8,4 = 15 


Fig. 1. Example of constructing binary tree based on sorting 


There is Theorem 1 [3]. For a single-dimensional array of N components, a binary tree can be built in parallel 
N’ 
using sorting with 7 (=) = O(log, N) time complexity. 


The used sorting is stable; the binary tree is implied to be constructed with uniqueness. The indices of all me- 
dial components (all roots of subtrees) can be identified [3]. Considering this modification, all the indices from the 
above example for N subtree values can be calculated synchronously and mutually independently. This leads to a sin- 
gle estimate of the build time of the binary tree. For each specific N , all the values of the tree node indices can be cal- 
culated a priori and stored in the computer memory. With their help, the sorted components can be synchronously and 
mutually independently addressed to all the tree nodes. Formulas for calculating the node indices depend only on the 
total number of N input elements and are in no way dependent on their mutual arrangement after the stable sorting. To 
simplify memory addressing, the computed indices can be ordered at each level and arranged in ascending levels. Then, 
the entire population of the ordered node indices is read from N key. It only remains to arrange the sorted tree elements 
by the read-in addresses. Based on the above, there is 


Theorem 2. For a single-dimensional array of N components, a binary tree can be built in parallel using sorting and 
N’ 
prior calculation of indices with T (x = O(1) time complexity. 


The following unified table contains the formal estimates of time complexity of sequential and parallel algo- 


rithms for constructing a binary tree versus the proposed algorithms. 
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Table 1 
Comparative estimates of time complexity of sequential and parallel algorithms for constructing binary tree versus 


proposed algorithms 


Acceleration when using 


Algorithm time Acceleration when using unit 


Binary tree algorithm logarithmic time-complexity 


complexity time-complexity algorithm algorithm 


Algorithm of A. 
Lagana and V. Kumar 
(2004) [5] 


Algorithm of P. 
Chalermsook (2015) 


[6] 
Polynomial algorithm 
(2016) [7] 


“Left child — right 
sibling” algorithm 
(2014) [8] 


Pattern-based 
algorithm (1991) [9] 


The presented 
algorithm with 
logarithmic estimate 
of time complexity 
(2015) [3] 

The presented 
algorithm with single 
estimate of time 
complexity (2015) [3] 





In Table 1: D is capacity of the template dictionary, N is the number of input elements of the binary tree, k 
is the dimension of the space in which sorting is performed. 


The table shows that the proposed algorithm with a logarithmic estimate of time complexity abstractly 1m- 
proves estimates of the known algorithms. Minimum acceleration is achieved with respect to the algorithm from [5]: 


~ 


S|) 


= O(N ° In2) ,or —= O(N te and maximum acceleration is achieved relative to the polynomial algorithm from [7]: 





T 

T N° ie : 

r =O ary or r = O(N ) . Regarding the proposed algorithm with a single estimate of time-complexity, the eval- 
Og, 


uation of the known algorithms also improves. In this case, minimum acceleration is achieved relative to the algorithm 


~ 


from [5]: r = O(N ) , and maximum acceleration is achieved with respect to the polynomial algorithm from [7]: 


* 


p= O(N’). 


Conclusion. The developed algorithms differ from the known techniques [5—7, 10, 11] of constructing a binary 
tree in that they use maximum parallel sorting to calculate the indices of the nodes. In this case, either a logarithmic 
number of steps is consumed by building a tree, or additional time is not spent at all, if the values of the indices are a 
priori calculated for all N values in some real boundaries and stored in the computer memory. The proposed parallel 
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algorithm for constructing a binary tree can be used to organize efficient methods for dynamic processing of databases. 
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